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TECHNICAL MEMORANDUM 


ANALYSIS OF LIGHTNING FIELD CHANGES PRODUCED 
BY FLORIDA THUNDERSTORMS 

CHAPTER 1 
INTRODUCTION 

Thunderclouds produce electric fields at the ground that can approach several kilovolts 
per meter, and are usually directed upward. Pointed objects near the ground distort and inten- 
sify these fields, and ensuing corona currents lead to the formation of a space charge shielding 
layer. Because of this shielding layer, ground observations of the electric field no longer 
contain enough information abcut the details of the thundercloud charges aloft. Fortunately, a 
great deal can still be learned about the electrical structure of thunderclouds by analyzing 
ground-based measurements of the changes in the electric field that are caused by lightning. 
Since the space charges cannot rearrange quickly enough during a lightning event (typically 0.3 
to 0.4 seconds), the effect of shielding can usually be ignored. 

C.T.R. Wilson undertook the first systematic study of lightning field changes [Wilson, 
1916]. He measured AE’s at various distances from thunderstorms and showed that distant 
cloud discharges produced field changes that were usually reversed in polarity from closer 
discharges. From these observations, Wilson [1920] deduced the classic bipolar model for 
thundercloud charges, i.e., positive charge at high altitudes and negative charge at lower 
altitudes [Wilson, 1929]. 

Unfortunately, Wilson’s measurements were limited to a single station and could not 
provide a detailed picture of the charges that were involved in lightning. In order to obtain 
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more information about how lightning altered the cloud charge distribution. Workman and 
Holzer [1942], Workman, Holzer, and Pelsor [1942], Reynolds and Neill [1955], and later 
Krehbiel [1979] made multiple station measurements of electric field changes in New Mexico. 
Krehbiel [1981] has recently reviewed the measurement techniques, the methods of data 
analysis, and the results. These authors found that many flashes could be described by simple 
charge models, and the resulting solutions were interesting not only because they provided 
quantitative estimates of the lightning charge transfers, but also because they provided an 
estimate of the altitudes of the charges (and corresponding temperatures) and the charge 
locations relative to radar echoes. The New Mexico results showed that lightning charges were 
located at temperatures that were below freezing, a result that has been verified by more recent 
work [Barnard, 1951; Hacking, 1954; Tamura, 1958; Hatakeyama, 1958; Takeuti, 1966; 
Ogawa and Brook, 1969]. This result in turn has motivated several important laboratory 
experiments on the role of ice particles in cloud electrification processes [Workman and 
Reynolds, 1949; Reynolds, Brook, and Gourley, 1957; and more recently Caranti and 
Illingworth, 1980; Jayaratne, Hallet, and Saunders, 1980; Latham, 1981; Illingworth, 1985; 
Baker et al., 1987; and Saunders and Zhang, 1987]. 

In the early 1970’s, the NASA Kennedy Space Center (KSC) installed a large network 
of electric field mills to identify atmospheric electrical hazards that might threaten the launches 
of spacecraft or ground operations. This network, located in the maritime tropical climate on 
the Florida peninsula, contained 21 ground-based sensors and covered a total area of 
approximately 15 x 25 km 2 . 
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Jacobson and Krider [1976] used the KSC network to study thunderstorm fields and 
developed a nonlinear least-squares minimization procedure (that will be described in detail 
below), to analyze lightning field changes. They found that discharges to ground effectively 
deposit positive charge (or, equivalently, remove negative charge) at altitudes of 6 to 9 km 
where corresponding ambient air temperatures are between -10 °C and -34 °C. Although these 
charge altitudes are somewhat higher than in New Mexico, the corresponding air temperatures 
were typical of other regions. Jacobson and Knder [1976] also found that a significant fraction 
of the discharges to ground produced field changes that were small or even reverse polarity 
within 3 km of the strike point. An analysis of these events showed that a small volume of 
positive charge, 0.4 to 5 C, was often neutralized at altitudes between 1 and 3 km, together 
with the main negative charge. 

Maier and Krider [1986] analyzed KSC field mill data taken during the summer months 
of 1976-1978 that was similar to the data used by Jacobson and Krider [1976]. These authors 
used a computer algorithm to identify lightning and compute values of AE, rather than manual 
methods used by Jacobson and Knder [1976]. The results indicated that the negative charge 
centers remained at constant altitude throughout the life cycle of the thunderstorm. The mean 
charge altitude depended somewhat on the storm, but the values ranged from 6.9 km (-14 C) 
to 8.8 km (-26 °C), with total charge transfers ranging from 11 °C to 44 °C. 

A few years later, Koshak and Krider [1989] developed an improved computer method 
for computing accurate field change values and extended the work of Maier and Krider [1986] 
to active thunderstorms. Koshak and Krider [1989] also analyzed intracloud lightning and 
showed that a typical cloud discharge produced moment changes of 1 13 °C km to 343 C km 
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and that high altitude cloud flashes were systematically larger than lower altitude flashes. 
Ground discharges in active storms were found to deposit charges at altitudes and ambient 
temperatures that were comparable to those in small storms, and the total charge values were 
also similar. Krider [1989] later pointed out that the apparent separation between the positive 
and negative charge centers depended on the storm size. 

One of the more interesting results of Koshak and Krider [1989] is the high degree of 
symmetry in the directions of cloud discharges. High altitude cloud discharges effectively 
transported positive charge in a downward direction toward the negative charge region, while 
low altitude discharges transported charge in an upward direction. Cloud discharges that 
occured at the same altitude as the negative charge region tended to transport positive charge in 
a horizontal direction. These results are consistent with the classic bipolar charge structure 
proposed by Wilson [ 1929] and with the extended tripolar structure that is reviewed in 
Williams [1989]. 

Although least-squares analyses of the KSC field mill data have greatly improved our 
knowledge of lightning charges, these studies have several limitations and inherent biases 
associated with them. First, the computer methods for detecting and computing lightning AE’s 
from field mill data are not perfect. Subsequent least-squares analyses of erroneous AE data 
can lead to nonphysical charge solutions. Second, since very small field change values have 
proportionately larger errors (see comments in Chapter 2 on digitization errors in KSC field 
data), only large lightning events have been analyzed, i.e., those that produce a AE > 1 kV/m at 
three or more field mill sites. Hence, the solution statistics are biased because the smaller 
events have not been included. Finally, there are biases because complex lightning events can 


4 



not be described with the simple charge models used. More complicated charge models have 
not been used because a finite set of measurements do not, by themselves, provide enough 
information to uniquely describe complex events. With a normal uncertainty in the AE values, 
there will always be an uncertainty in the charge solutions. Even without measurement error, 
the nature of Coulomb’s Law is such that there will always be an infinite number of charge 
distributions that will produce exactly the same set of field values. In order to determine a 
"correct" solution to an inherently non-unique problem, we will need to add extra physical 
constraints so that the number of our choices are reduced. Such constraints should be based on 
our knowledge about the physics of lightning and the thunderstorm environment. 

In the following chapters, we will describe an interactive computer program that 
combines manual and automatic analysis techniques to improve lightning identification and the 
computation of AE values (regardless of their amplitude). This new program has been applied 
to eight Florida thunderstorms and the resulting sets of AE values have been analyzed using the 
standard least-squares approach and point charge (Q) and point dipole (P) models. To reduce 
biases from a AE threshold requirement, we have analyzed all flashes that produced a field 
change of 1 kV/m at two or more sites (rather than three or more sites as described above). 
Since we have analyzed the entire life cycle of each storm (some lasting 2-3 hours), and since 
the interactive program has allowed us to compute improved AE values and has given us 
confidence in reducing our AE-threshold requirement, we have obtained improved Q- and 
P-statistics over previous studies. 

The inability to describe complex lightning events with simple charge models is a far 
more fundamental problem. Because of this, we have also developed an entirely different 
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approach for analyzing lightning field changes. This new method can be used to describe very 
complicated lightning charge distributions and provides an easier way of adding arbitrary 
external constraints to the solution process. Solutions are no longer the parameters of simple 
charge models, but are volume distributions of charge that are defined on a grid of finite 
dimension and resolution. In order to describe lightning events having a variety of locations, 
the grid is centered above the KSC network and has a large detection volume. In addition, our 
method provides a framework where we can use a standard eigenanalysis to determine the 
information content of the field change measurements, and assess the effects of measurement 
errors, network geometry, and solution grid geometry on the accuracy of the solution. 

The chapters below are organized in the following way. All details of the KSC 
measuring network, examples of thunderstorm fields, and a discussion of the interactive AE 
program are provided in Chapter 2. In Chapter 3 we review the least-squares approach of 
Jacobson and Krider [1976] for analyzing AE’s, and describe in greater detail some limitations 
in the method. We next introduce our new analysis method by showing that we can describe 
field changes in terms of a linear system of equations. Various algorithms for solving this 
system will be discussed near the end of Chapter 3, and we will note that the most appropriate 
algorithm is the "method of steepest descent". This algorithm allows arbitrary external 
constraints to be added to the solution process and the importance of these constraints will be 
stressed in Chapter 4, where we give several examples. 

As a check of our new method, we analyze the fields produced by known lightning 
sources in Chapter 5 and discuss solution errors. We show that a particular form of the method 
of steepest descent, a Landweber iteration, is most useful in obtaining accurate solutions. We 
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give a theoretical description of solution errors and we also explore the information content of 
the data using an eigenanalysis approach. In Chapter 6, we give all of our Q- and P-model 
results for eight Florida storms. We also analyze three lightning events from one of these 
storms using the Landweber iterative technique and then compare the results with the 
associated Q- and P-model results. Chapter 6 ends with some suggestions for future work, and 
Chapter 7 summarizes our results. 
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CHAPTER 2 


MEASUREMENTS 

This chapter describes our measurements of thunderstorm electric fields and a computer 
algorithm that has been used to determine the values of lightning-caused changes in these 
fields. The finite set of derived field changes will loosely be referred to as "measurements" in 
later chapters. The net errors in the field change values are important in determining the 
accuracy of our solutions (see Chapter 5 on solution error). Accordingly, this chapter also 
reviews various difficulties associated with computing accurate AE’s from the electric field 
data and concludes by giving reasonable estimates of our overall error. 

2.1 Electric field measurements 

The electric field data have been obtained at the NASA Kennedy Space Center (KSC) 
and the Cape Canaveral Air Force Station (CCAFS). These facilities are located on the 
east-central coast of the Florida peninsula. During the summer months, thermally driven 
sea-breeze circulations develop in this area, interact with the Westerlies, and produce lines of 
convection that often produce heavy precipitation and lightnning. In order to identify atmos- 
pheric electrical hazards, KSC and CCAFS operate and maintain a large ground-based network 
of electric field mills at the sites shown in Figure 2.1. 

Figure 2.2 shows a photograph of a typical field mill site. Note that each sensor is 
placed on a level, paved surface that is free from grass and other protruding obstacles, such as 
trees, buildings, etc., that might distort the local field. The sensors are cleaned at monthly 
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intervals so any adverse effects from the elements (e.g., sea spray salt deposits, insects, etc.) 
are minimized. 

Each field mill measures the vertical component of the atmospheric electric field, and 
during normal operation, the data are sent to a central recording station, where they are 
digitized at a rate of 10 samples per second and stored on nine-track magnetic tape. The 
digitization accuracy is about 30 V/m and each field mill has a dynamic range of -15 kV/m to 
+15 kV/m. In order to remove 60 cycle fields, the sensors operate with a low pass filter that has 
a time constant of about 0.1 seconds. Although this constant is too slow to time-resolve the 
individual components of a multiple stroke flash, it is more than adequate to resolve the entire 
discharge (a typical discharge has a duration of about 0.5 seconds). Absolute calibration of 
each sensor has been estimated to be accurate to within about 10% [Jacobson and Krider, 

1976]. All of the field mill data that are analyzed in this study were obtained in 1978 as part of 
the Florida Thunderstorm Research International Project (TRIP) [Pierce, 1976]. 

2.2 Computation of lightning-caused field changes 

Computer methods for identifying lightning and computing values of lightning-caused 
changes in the electric field have been described by Maier and Krider [1986] and Koshak and 
Krider [1989]. In this section, we will review these methods and discuss their deficiences. 

With these deficiencies in mind, we will then introduce an improved interactive AE program. 
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2.2.1 Identifying lightning in field mill data 

Figure 2.3 is an example of the electric field records that were obtained simultaneously 
at three different field mill sites. These data were obtained at the NASA Kennedy Space 
Center during a fairly active thunderstorm on July 6, 1978. In this figure, a positive potential 
gradient is the polarity that would be produced by a positive charge aloft, or equivalently, a 
negative charge aloft will produce an upward directed field. The abrupt field discontinuities 
throughout the record are produced by lightning flashes, and it is these field changes that need 
to be identified and computed. 

Since individual Florida thunderstorms can produce thousands of discharges 
[Livingston and Krider, 1978], manual detection and analysis of each individual lightning 
event is not a practical procedure. The most obvious characteristic that distinguishes a 
lightning event from other variations in the field is the simultaneous occurrence of a large field 
derivative of either plus or minus polarity that lasts a few tenths of a second at two or more 
field mill sites. Accordingly, the computer algorithms developed by Jacobson and Krider 
[1976], Maier and Krider [1986], and Koshak and Krider [1989] have used this characteristic 
for lightning identification. Jacobson and Krider [1976] used an algorithm that began by 
arranging the field mill data into consecutive 1 -second time blocks. The averages of the first 
two field values in each adjacent block were used to extrapolate a field value to the end of the 
third time block. If the actual field at the end of the third time block differed by more than 600 
V/m from the extrapolated field at two or more sites, a lightning was assumed to have 
occurred. The time of the flash was taken to be the time at which the majority of sites 
experienced their maximum field derivative. 
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2.2.2 Obtaining initial and final fields 

Once a lightning event has been found, the time the discharge began (h) and the time 
that i, ended W need to be determined in order to compute a fieid change. Examples of these 

times are shown in Figure 2.3 for a large flash that occurred near field mill site 2 just pnor to 
190700 GMT. Maier and Krider [1986] found the beginning and ending times by calculating 
the field derivative before and after the flash time and then assigning h and t, to be the times 

when most sites had field derivatives that wete consistently below 1750 V/nVs. Because the 
shapes of the field changes are often very different from site to site for a particular lightning 
event, Koshak and Krider 119891 improved this procedure by using a a pattern recognition 
algorithm to move time markers to the proper times tj and t, . 

2.2.3 The time-varying background field 

Since the background electric field is often changing when a discharge occurs, one 

cannot simply subtract the initial field E(t,) from the final field E(« r ) to obtain an accurate value 

of 4E. Note in Figure 2.3 that the flash just before 190700 GMT at field mill site 2, is 
preceded by a large field derivative just before t, and just after h. and a similar behavior can be 

seen for many other flashes in Figure 2.3. These field variations are caused by currents within 
and near the thundercloud, and the associated displacement current densities (i.e„ e„dE/dt) can 

approach tens of nanoamperes per square meter [Krider and Blakeslee, 1985; Koshak and 
Krider, 1989; Denver and Krider, 1990], More accurate field change values can be found 

using: 
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A E= E f -^\ f (Aj/2) - £+^|.(Af/2) 


( 2 . 1 ) 


In this procedure, about six field values before tj and after tf were used to calculate the field 
derivatives, and At = tf - tj was the duration of the lightning event. Koshak and Krider [1989] 
have shown that the correction term, (- dE/dtl f - dE/dtl j) At/2, can be as large as 40% of the 

uncorrected field change during active thunderstorms. 

The pattern recognition algorithm used by Koshak and Krider [1989] helped to 
determine a more accurate tj and t f for each flash. This procedure also provided more accurate 

field projections. If for example, tj and tf were only approximated, then the interval. At, and the 
field values used to calculate dE/dtl j and dE/dtl f could be in error. In a worst case scenario, if 
two lightning events are very close to one another in time, a slight error in tj and tf could 
produce a background correction that, in reality, is associated with another lightning event. 

2.2.4 L- and F-changes 

Another factor that needs to be considered in the computation of field changes is the 
possible presence of a leader L- (leader) or F- (final) change structure in the shape of AE 
[Pierce, 1955]. In Figure 2.3, site 2 contains a F-change structure at about 190516 GMT. In 
order to show this more clearly, Figure 2.4 gives an expanded view of the flash of interest. 

Note the initial sharp discontinuity that defines the start of the flash, and then the large positive 
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peak. After the peak there is a substantial negative field transition (or F-change) lasting about 
0.2 seconds. After the F-change, the field returns to a slowly varying background level. Since, 
we are interested only in the total change in the background level, we must ignore any L- or F- 
changes in AE. To detect these changes, Koshak and Knder [1989] check for large changes in 
the field derivative before and after the main field discontinuity. Again, the success of this 
approach depends on the success of pattern recognition iterations in finding the times directly 
before and after the main field discontinuity defining the flash. 

2.2.5 Detection of poor field mill data 

Unfortunately, the algorithm discussed by Koshak and Krider [1989] cannot be used to 
analyze poor field mill data. Examples of poor data are given in Figure 2.5, an expanded 
portion of the data given previously in Figure 2.3. Note how the output from mill 25 is being 
affected by a 1 Hz background noise signal, although it does still detect lightning events. Mill 
18 has spurious field variations that could be due to a local rain shower. In both cases, the 
lightning field variations, such as is recorded "cleanly" at site 2, is being distorted. 

At least two problems are caused by noise in the field mill data. The first is that the 
noise can satisfy the lightning detection criteria and produce a false event. As discussed above, 
a lightning signature is distinguished by a large field change at two or more field mill sites; 
thus two noisy sites can create a false event. The second problem arises when noisy data are 
used to compute the time varying background field. Koshak and Krider [1989] have 
minimized this problem by extrapolating only those fields that exibit a high degree of linearity, 
but this method is not perfect. 
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The records in Figure 2.3 show that data are missing at about 190930 GMT (see the 
long straight lines connecting the field values) probably due to a short power interruption. 
Automatic analyses of these records could result in erroneous AE values. 

2.2.6 Interactive algorithm 

In order to improve the detection and AE computation of lightning flashes and to 
eliminate bad field data, we have developed an interactive computer algorithm. The procedure 
starts by plotting the data on a high resolution CRT display and then making a visual 
examination to eliminate all poor field mill data. Typically, the computer displays about 10 
minutes of data for five field mill sites, with a temporal resolution of 10 samples per second. 
When searching for bad data, the operator can examine the next 10 minutes of data by simply 
pushing a mouse key. By scrolling through the data in this fashion, most noise and data drop 
outs can be detected and then eliminated from further analysis. 

After the bad data are eliminated, the operator returns to the beginning of the data file 
and then scrolls through the good data in order to identify individual lightning events. To do 
this, the operator selects three or more field mill sites that are spatially separated and then 
scrolls through the records and marks the times of sharp discontinuities that define each flash 
using the mouse. The amplitude of the field is automatically normalized on each data frame. 
After the operator identifies a flash time, the algorithm automatically calculates the values of 
AE using the pattern recognition iterations, L- and F-change checks, and the time varying field 
corrections that were described above. Within a few seconds, the initial and final times of the 
flash are plotted on the screen together with dotted lines that show the time-varying field slope 
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corrections. After the field changes have been computed, the flash times and AE values are 
stored on the computer hard disk. 

With this interactive program, the operator can make a decision as to how well a 
particular AE has been analyzed or if in fact the event should be analyzed manually. Figure 2.5 
shows a flash at about 190745 GMT (mill 2) that has a very complicated F-change structure. 

In this case, the interactive algorithm would not be able to obtain an accurate AE value with the 
normal F-change checks, so a manual analysis would be used. The algorithm helps the 
operator to carry out the AE computations, stores the values, and other housekeeping activities 

that would otherwise be very time-consuming. 


2.3 Errors 

Since ultimately we will use the AE values to compute the locations and magnitudes of 
changes in the cloud charge distribution, it is important to determine what effect errors in AE 
have on our final charge solution. Intuitively, we might expect that large uncertainties in the 
AE values will produce large uncertainties in the inferred solution, and as shown in section 5.1, 
this is a correct statement. What is less obvious, however, is that solution error does not 
necessarily increase in a one-to-one sense with measurement error. Depending on the location 
of field sensors for instance, even very small amounts of error in AE can result in huge solution 
errors. The effects of AE error on solution error are discussed more thoroughly in Chapters 4 

and 5. 

For the current discussion, we only need an estimate of the total error in AE. In our 
particular problem, the author feels confident that a reasonable upper bound to the AE error is 
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about 10% or 30 V/m, whichever is larger. This estimation includes all sources, e.g., random 
error in the field measurements, digitization error, inaccuracies in time-varying field and 
pattern recognition calculations, etc.. 
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CHAPTER 3 


FIELD CHANGE ANALYSES 

In this chapter, we will review the least-squares optimization procedure that has 
previously been used to analyze lightning field changes, and we will point out some basic 
limitations in the method. We will also introduce a new, linear method that can provide new 
opportunities and new insights. 

The primary aim of both analyses is the same. When a lightning discharge occurs, the 
original charge distribution in and around the cloud is altered, and we would like to describe 
the location and magnitude of this change. According to Coulomb’s Law, any change in the 
thundercloud charge will produce a change in the field at the ground. In textbook problems, 
one is usually asked to find the field that is produced by a given charge distribution in the 
presence of various conductors and/or dielectrics. Here we are interested in the inverse 
problem, i.e., we want to find the changes in the cloud charge distribution when we are given 
the values of AE. We have previously pointed out that this inverse problem is fundamentally 
non-unique. In Chapter 4, we will show how external constraints can be used to reduce the 
solution ambiguities. 


3.1 Nonlinear least-squares optimization 

Least-squares optimization methods were first applied to the analysis of lightning field 
changes by Jacobson and Krider [1976]. The analysis begins by assuming that the measured 
field changes, AE i; can all be described by a simple charge model. The parameters of the 

model (e.g., charge location and magnitude) are then inferred from the measurements by 
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minimizing a (reduced) Chi-square function of the form: 


, 2 = iyf AE, -mA 2 _ C 2 
r ^ c, J N f ' 


(3.1) 


where Mj = M^a*), k=l,..,p, is the model field change value at the i* field mill site; there are p 
model parameters a k ; a ; is the rms measurement error, and N f = m - p is the number of degrees 

of freedom (i.e., number of measurements minus the number of unknown model parameters). 
In practice, the model M, is nonlinear so nonlinear methods are used to search for the model 

parameters that minimize C r 2 . 

Because of measurement errors, the minimum values of C 2 for different events are not 
unique numbers, but are a set of numbers that are distributed according to the C 2 distribution. 
The expectation value of this distribution is simply the number of degrees of freedom in the 
solution (i.e., <C 2 > = N f ). When N f is increased, the standard deviation of the C 2 distribution 

decreases so that large values of C 2 become more improbable [Bevington, 1969]. In general, 
Jacobson and Krider [1976], Maier and Krider [1986], and Koshak and Krider [1989] have 
found that a value of C r 2 < 10 usually corresponds to a suitable solution when = 0.05. 


3.1.1 The Marquardt algorithm 

The nonlinear search for the optimum parameters has usually been an iterative 
algorithm first proposed by Marquardt [1963] and subsequently described by Bevington 


[1969]. In this algorithm, the model function Mj is linearized by expanding it in a Taylor series 
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and dropping all higher order terms. If we insert the linearized model function back into Eq. 
(3.1), C r 2 can now be regarded a function of the model parameter increments, Aa k : 


Cf(Aa k ) = 


^E— 2 (a£/- Mio 
N, M of V L 


p 

+ E 


<*Mj isa k 
da k o 


(3.2) 


In order to find the optimum values of the increments, i.e. the values that minimize C r 2 , we take 
derivatives of C r 2 with respect to each Aa k and set the results equal to zero. This produces p 
equations with p unknowns (parameter increments) that are of the form: 


E— 2 (A£, - M i0 ) 
■=i o, 


dMj 

da k 


= E A*, Afli 


1=1 


(3.3) 


where A u are the elements of an approximative curvature matrix (pxp) given by: 


V 1 dMi | dM , i 
A ki=^~l 


i = l of dtf* 1° 5^1 


(3.4) 


Thus, we can start with an initial guess of the model parameters and then iterate to find better 
values by solving Eq. (3.3). The procedure is repeated until a suitably small value of C r 2 is 

obtained. 

Solving Eq. (3.3) for the parameter increments is known as the expansion method; and 
this approach is particularly efficient when we are close to a minimum in the C r 2 function, but 

slow when we are far from a minimum. The Marquardt algorithm overcomes this deficiency 
by replacing the matrix A' with A = A' + XI, where X is a positive scalar and I is the identity 
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matrix. Now, the iteration formula can be written as: 


Aa = (A' + M) » V, ( 3 . 5 ) 

where V is a k-dimensional vector whose components are given by the terms that are summed 
on the left side of Eq. (3.3). At the start of the search, X is made large so that the diagonal 
elements of A dominate. The parameter increments in this case are in a direction nearly down 
the gradient of the C r 2 hypersurface and the minimum is approached quite rapidly. Since the 

gradient tends toward zero, and frequently changes direction as a minimum is approached, the 
value of X is decreased to favor the expansion method. 


3.1.2 The Q- and P- models 

The simplest charge model and one that can be used to describe a spherically symmetric 
lightning event is the point charge or simply the "Q- model" [Jacobson and Krider, 1976; Maier 
and Krider, 1986]. This model contains four unknown parameters, the charge location (X,Y,H) 
and magnitude AQ. With the Q-model, the field change at the i* field mill site is: 


M . - 2AQH 

j t 

4neo(H 2 + £>?) 2 


(3.6) 


where D 4 2 - (X r X) 2 + (Y r Y) 2 and (X^Yj) is the location of the i 1 * 1 field mill site. This model 
assumes that the field sensors are located on flat, perfectly conducting ground. The factor of 2 
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in the numerator of Eq. (3.6) comes from an image charge that satisfies the boundary condition 
at the ground. 

A model that is often used to describe a cloud discharge [Koshak and Knder, 1989; 
Krider, 1989] is a six-parameter point dipole or "P-model." Here, the discharge is assumed to 
produce field changes on the ground that are described by. 


Mi=- 


1 

2k 


f 3//(APR,) 



(3.7) 


where the point dipole vector AP has the components (AP X , AP y , APJ, and is located at 
(X,Y,H). The position vector, R ; , points from the i* field mill site to the point dipole. Figure 

3.1 shows the geometrical aspects of both the Q- and P-models. In practice, each model is run 
on the same lightning event and the model that produces the lowest value of C r is 

assumed to be correct (provided, of course, that C r 2 < 10) [Maier and Knder, 1986]. 


3.1.3 Limitations of model fits 

A large fraction of lightning discharges alter the cloud charge in ways that produce 
complex field change patterns, e.g., multi-branched air discharges. Under these circumstances, 
it is usually not possible to describe the event using a simple Q- or P-model (Koshak and 
Krider [1989] could describe only about 50% of all large lightning events using Q- and P- 
models). A reasonable response to such a problem might be to invent more complicated charge 
models that have several model parameters. As will be seen below, and again in Chapter 4, 
this is not an acceptable course of action. 
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In the studies done by Workman, Holzer, and Pelsor [1942] and Krehbiel [1981], it was 
suggested that the number of measurements sets an upper limit to the complexity of the charge 
model used. In effect, it was asserted that a solution could not be found because the number of 
measurements (m) were fewer than the number of model parameters (p). 

On the contrary this is not the reason for being restricted to simple models. For 
instance, one can find the minimum (^(a) when m < p by simply checking all reasonable 
values of the parameters until the minimum is found. This can be accomplished on the 
computer by using nested DO-LOOPS with adequate resolution for parameter increments. 
Furthermore, the reader is urged to note that the Marquardt algorithm, as it stands, can also be 
used to find the minimum when m < p. Note that for m < p, the approximative curvature 
matrix A given above is singular (i.e., noninvertible). However, as pointed out by Marquardt 
[1963], the addition of the term XI to A removes this singularity so that the minimum can be 
approached (i.e., A = (A' + XI) is invertible). Bevington [1969] does not mention this point, 
and in fact, Marquardt [1963] devoted only one sentence to it in his original discussion (see 
Marquardt [1970] for a full discussion). Thus, even when m < p the measurements allow one 
to find a solution; this solution is, as always, non-unique. 

With these comments in mind, the primary reason for avoiding many parameter models 
is due to the difficulties associated with adequately constraining the parameters so that physi- 
cally reasonable solutions are obtained. Unfortunately, the measurements alone do not always 
provide enough constraint to the model parameters so that more than one solution may be 
deemed acceptable. In effect, there exist multiple minima in the C 2 hypersurface; if additional 
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(or "external") constraints are not added to the problem, no clear cut determination of the 
"proper" solution can be made. 

In all least-squares analyses to date (regardless of the algorithm used to find the 
minimum C r 2 ), few constraints, if any have been applied to model parameters. The Marquardt 

algorithm, most often used, is no exception and this is clearly demonstrated in section 4.2 
below. Note that there is the possibility of redefining the error function to be C r 2 [given in Eq. 

(3.1)] plus some additional terms that directly constrain the model parameters, however a more 

desirable approach will be described in section 3.2. 

Another limitation of the least-squares methods used to date is that there is no practical 
way to determine the information content of the measurements. We will now introduce a new 
analysis method that allows a standard eigenanalysis to be used to determine information 
content and which provides other advantages as well. 

3.2 Fredholm integral formulation 

In order to overcome some of the limitations of the least-squares model approach, we 
will reformulate our inversion problem in terms of a linear system of equations. This is a new 
and fundamentally different approach to the analysis of lightning field changes. 

Many inversion problems have previously been put into the form of a linear Fredholm 
integral equation of the first kind [Twomey, 1977]: 

g(u) = \K(u,v)f(v)dv , ( 3 - { 
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where g(u) represents measurements, f(v) is an unknown distribution that is to be inferred, and 
K(u,v) is a kernel function relating the measurements to the unknown distribution. Several 
examples of physical problems having this general form are given by Twomey [1977]. For 
instance, an atmospheric temperature profile T(z) [= f(v)] can be inferred from measurements 
of the radiation intensity I(X) [which, apart from a few constants, is g(u)] at different 
wavelengths X [= u], the kernel being a derivative of an optical transmissivity function. In this 
case Eq. (3.8) becomes a modified form of the equation of radiative transfer. We will soon see 
that Coulomb’s Law can also be written in the form of Eq. (3.8). 

The formal redevelopment of the problem begins with an integral description of the 
electrostatic potential due to a known volume distribution of charge (see Jackson [1975]). 
Starting with the well-known divergence theorem: 


j V' AdV' = j Anda' (3.9) 

V' s' 

if we let A = <|>V'G - GV'<J>, where <{> and G are arbitrary functions of space, we get: 

J ($V' 2 G-GV' 2 $)dV' = j($V'G-GV'ty) h da' (3.10) 

V' s' 

which is Green’s second identity or Green’s theorem. If we identify <J> and G as the scalar 
potential and a Green function, respectively, we have the following constraints: 
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V' 2 <t )(r') = -p (r')le 0 


(3.11) 


V' 2 G(r,r') = -8(r-r')/% 


(3.12) 


G(r,r ') 


1 

4;teo |r-r'| 


+ F(r,r'). 


(3.13) 


Here, the Green function is simply the potential of a unit point charge at r' plus the potential 
F(r,r') due to a system of image charges outside V'. Substituting Eqs. (3.1 1) and (3.12) into 
(3.10) and carrying out the volume integration over the delta function, we obtain an integral 
relation for the potential: 

(J>(r)= j G(r f r')p(r')dV'-eo (3.14) 

To solve Eq. (3.14) for the potential, we can pick F(r,r') so that either Dirichlet [G(r,r') 
= 0 on S'] or Neumann [3G(r,r')/5n' = 0 on S'] boundary conditions are satisfied depending on 
whether we know the value of <(>(r') or dfy/dn' on S', respectively. For our purposes, the upper 
half space is the volume of interest, V', and we require that the potential be zero everywhere 
on S' as shown in Figure 3.2. In this Dirichlet problem, the integral relation for the potential 
reduces to: 

<Kr)= J G(r, r')p(r') dV'. (3.15) 

V' 
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Using the method of images, the appropriate form of F(r,r') can be found so that the 


boundary conditions are satisfied. The potential in V' due to a unit charge in V' and its image 
below the conducting plane is: 

C(r ’ r ' >= 4i5(Vr} < 316 > 

where R = [(x-x') 2 +(y-y') 2 +(z-z') 2 ] 1/2 and T = [(x-x') 2 +(y-y') 2 + (z+z') 2 ] 1/2 . Substituting (3.16) 
into (3.15) gives: 

(3.i7) 


At this point, we can take the gradient of Eq. (3.17) to determine the electric field at all points 
within V'; if a lightning discharge alters p(r'), we can calculate the temporal change in the 
field (or potential gradient) at the ground using: 


A(Vf(r) ^ ) = 


z f 2z" Ap(r') 

4 *eo l ft*-*') 2 +(y-y') 2 + z^] 3/2 


dV'f 


(3.18) 


where Ap(r') describes the change in the cloud charge distribution that was produced by the 
lightning discharge. 

Adopting the notation of Twomey [1977] we may denote: 


g(x,y) = A(V<[>(r) ^ -z 


(3.19) 
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K(x,y,r') = z'/(2ke 0 l(x-x'f + (y-y'P + z' 2 ] 312 ) 


(3.20) 


f(r') = Ap(r') . (3.21) 

With these changes in the variables, equation (3.18) becomes: 

g(x,y) = J K(x,y,r')f(r')dV' (3.22) 

which is essentially the same as the Fredholm integral equation given previously in Eq. (3.8), 
with (x,y) — » u and r' — » v. Here, g(x,y) describes the changes in field anywhere on the 
ground, f(r') is an unknown charge density distribution, and K(x,y,r') is a geometrical kernel 
function relating the two. Given a particular (x,y) field mill location on the ground, the 
magnitude of the kernel function K(x,y,r') can be plotted for arbitrary points in space r'. 
Figures 3.3 and 3.4 show such plots for representative distances D = [(x-x') 2 +(y- y') 2 ] m from, 
and altitudes z above, an arbitrary field mill site location (x,y). 

Before proceeding any further, it will be worthwhile to express Eq. (3.22) in matrix 
form. We can write the unknown source function as a series of (arbitrary) discrete charges at 
known locations on a grid as 

n 

fir ') = £ fj 5 (r ' ' - r ; ) , (3.23) 

7=1 
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where rj is the position of the j A grid point and fj is the value of the change in charge at that 
point. The field changes at position (x,y) on the ground can then be described by: 

ft 

80c,y) = 'LK(x,ys j )fj . (3.24) 

Since our measurements are made at a finite number of discrete points on the x-y plane, we can 
also write (3.24) as: 

n 

gi =U Kijfj i = 1 ,../n , (3.25) 

where m is the number of field mill sites, and n is the number of grid points in the upper half 
space. The matrix form of Eq. (3.25) is simply: 


g = Kf , (3.26) 

where g is now a (mxl) column vector of m field change measurements, f is a (nxl) column 
vector of the changes in charge at n grid points in the upper half space, and K is a (mxn) kernel 
matrix relating fields to charges. 

Thus, by using only Coulomb’s Law and the divergence theorem, we are now able to 
describe our measurements in terms of a linear system of equations. The nonlinearity is 
present only in the kernel functions as seen in Figures (3.3) and (3.4). Note that these 
functions are determined by simply evaluating the gradient of the Green function at the ground. 
Equation (3.26) is completely analogous to the temperature inversion problem that we 
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mentioned in section 3.1.3 and should not be confused with the linear system given previously 
in Eq. (3.3). 

A number of methods have been developed to solve linear systems such as Eq. (3.26), 
and many have been discussed by Twomey [1977]. Depending on the particular physical 
problem at hand, some of the methods of solution are undoubtedly more appropriate than 
others. 

An important and unavoidable characteristic of our approach is the fact that we now 
have a large number of grid points in the upper half space to consider and we must find a 
charge value at each of these points. In the temperature inversion discussed above, T(z) is 
described by using a single column of grid points, but we have the problem of finding f(r') = 
f(x',y',z') on a cubical (or cylindrical, etc.) volume of grid points. Since the number of 
columns in the kernel matrix K is equal to the total number of grid points in the volume, we 
can expect the size of the kernel matrix to be rather large. Methods of solving Eq. (3.26) that 
avoid the need to store and invert a large K matrix are described below. 

3.2.1 Gradient-constrained linear inversion algorithm 

A standard procedure for solving a matrix equation like Eq. (3.26) is the method of 
constrained linear inversion [Twomey, 1977]. In this procedure, one starts with an error 
function of the form 

e(f) = (g-Kf) 2 + jfHf, (3.27) 
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where the tilde symbol denotes transposition, H is a constraint matrix, and y is a factor that 
weights the constraint term /Hf. If we now take the derivative of e with respect to each of the 

unknown charge values, fj (j = l,...,n) and set the result equal to zero, we obtain n equations in 
n unknowns that have the solution: 

/ = (KK + yH ) -1 Kg . (3.28) 

Basically, this inversion formula provides values of f that minimize the error function, e. 

As noted above, we would like to have a procedure that avoids the need for storing or 
inverting a large kernel matrix. In Eq. (3.28), the matrix to be inverted is intolerably large, 
because it has dimension (n x n), where n is the total number of points in our volume grid. One 
alternative is to select only a few source points (i.e., limit n to 1 or 2 or 3, etc.) and then use an 
iterative procedure to move the sources around until an acceptable value of e is obtained. This 
procedure is conceptually similar to the nonlinear search procedure discussed in section 3.1.2, 
except that now we have an error function with constraints on the charge parameters, f. 

In order to find an optimum charge solution we will introduce an n-dimensional vector, 

P = (x 1 ,yi,z 1 ,...,x n ,y n ,z n ), that describes the locations of n grid points above the 

conducting plane with respect to an arbitrary origin. An iterative procedure that can reduce the 
error function to acceptable levels is to use the following steps: 

(a) choose a starting p, 

(b) find the optimum f using Eq. (3.28), 
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(c) form a new p' = p - tVe (f fixed), 

(d) return to step (b), 


(3.29) 


where the gradient operator V = (d/dx^d/d^), and t is an adjustable parameter. In practice, 
we have obtained satisfactory results using this procedure by choosing t = lApl/IVel = (Ar^ 

+...+ Ar n 2 ) 1/2 /lVel = (nd 2 ) 1/2 /IAel, where Arj 2 = Axj 2 + Ay, 2 + Azj 2 and d = 100 meters. 

Basically, one starts with an initial guess of the grid point locations, finds f by matrix 
inversion, and then further decreases the error function by moving the grid points down the 
gradient of e (holding f fixed). When the new grid point locations are found, the process is 
repeated. After a few hundred iterations, a much improved value of e is obtained. Since this 
process is fairly fast, and since the hypersurface e(p) may contain many local minima, we 
normally choose 1000 or more initial grid point configurations at random, followed by steps 
(b) through (d) until a suitable value of e is obtained. If all these randomized grid point 
configurations do not produce a sufficiently small value of e, the value of n is increased by one 
and the whole procedure is repeated. 

Note that because the term y H will dominate small eigenvalues in K K (see section 

4.1), the value of n can, in principle, be larger than the number of measurements. In addition, 
since e and f are only functions of p, the constraints we have imposed on f (through H and n) 
reduce the entire problem to one of finding the optimum vector p. 

In principle, the above iteration could be done starting with any initial charge 
configuration. Intuitively, however, one might expect to obtain better results if the first guess 
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places the source points fairly close to the lightning discharge. To acheive this, we first 
assume that the discharge is localized about some point r above the ground. Equation (3.26) 
then becomes: g = Kf = f[K,(r]) + ... + ^K^r,,), where the vectors K,^) are the column 

vectors of the kernel matrix. Since the flash is localized, r,= ... = r n = r, so that Kj(r,) = — = 
K„(r n ) = K(r). With this approximation, 

g = (f 1 + -+f n )K(r) = QK(r) (3.30) 

which is identical to the classic Q-model result given in Eq. (3.6). Now, to obtain an optimum 
value of Q(r) that is not too large, we can construct an error function of the form: e = [g - 
QK(r)] 2 + y Q 2 which is a degenerate form of Eq. (3.27) with H = I (identity matrix) and n = 1. 
The optimum value of Q at r is then given, in analogy to the constrained inversion result given 
in Eq. (3.28) as: Q opt = E[K i (r)g i ]/[E i (K i 2 (r)) + y]. Substituting this into our expression for e 

gives the minimum error e opt (r) at any point above the ground. By scanning the upper half 
space on a grid with 2 km resolution, we can find the value of r for which e^r) is a minimum, 
and then this value of r (denoted by r min ) becomes an optimum starting point (in the least 

squares sense). In practice, we can obtain a final solution by selecting initial charge 
configurations that are confined to a cube whose dimensions are 10 km on a side, and whose 
center point is at r,,^. The advantages of the gradient constrained linear inversion over the 

original Marquardt procedure are discussed in Chapter 4. 
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3.2.2 Twomey-Chahine algorithm 


This method, originally proposed by Chahine [1970] and later modified by Twomey 
[1975], is an iterative procedure that can be used to solve a matrix equation like (3.26) when 
the number of unknowns is large. As in any iterative approach, we start with an initial guess 
of the unknown distribution f°, and this guess is thereafter improved until the residual (g - Kf ) 2 
is sufficiently small. In tomographic approaches to solving two-dimensional remote sensing 
problems, Twomey [1987] reports successful reconstructions of up to 40,000 unknown 
quantities. The iteration scheme used to solve these highly underdetermined systems is given 
by: 

fj mX =/f(l + Cf^)> (3-31) 

where K- 5 = K ij /(K ij ) max is a scaled kernel element and Ci N = (g / J K^Of^rldV' - 1), 
though other forms of are possible. Twomey [1975] has shown that successive iterations 

by Eq. (3.31) lead to a solution that is constrained to the space spanned by the kernel functions, 
i.e., f is a linear combination of the kernel functions in the form f(r') = Zo^K^r'). Note that all 

portions of f(r') that are orthogonal to Kj(r') are not "seen" by the measurements and therefore 

are not constrained by the measurements. Since f(r') is forced to be nonorthogonal to the 
kernel functions, the solution is usually sufficiently constrained by the measurements so that a 
stable (nonoscillatory) solution results. We have verified that the results are stable when the 
source is a single point charge lying above the field mill network. 
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3.2.3 Method of steepest descent 


Unfortunately, the Twomey-Chahine iteration was designed to retrieve unknown 
distributions that are only one polarity (usually positive). We do not need this "polarity 
constraint" in our problem, since lightning can obviously involve both positive and negative 
charges. In addition, since the solution is confined to the space spanned by the kernel 
functions, it is difficult to add other contraints to the solution. In order to overcome these 
difficulties, we have devised a more general iterative procedure that allows the unknown 
distribution to take on positive, negative, and zero values, and that also allows for arbitrary 
external constraints. 

We define an error function to be of the form: 

e(f) = (g - Kf? +jjCj(f) + y 2 C 2 (f) + ... , (3.32) 

where 0,(0, C 2 (f)... are the constraint terms, and y x , y 2 ,... are the associated weighting factors. 
The iteration then becomes: 

f' = f-tVe, (3.33) 

where t is an adjustable parameter that determines how far one steps down the gradient of the 
hypersurface e(f), and where the gradient operator V = ? , (9/3fj) + ... +f n (3/3f n ). In essence, 

Eq. (3.33) is nothing more than a method of steepest descent applied to a constrained 
least-squares problem. Fig. 3.5 illustrates the basic idea behind this method. For clarity, the 
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illustration is for the simple case n = 2, i.e. the column vector f has only two components. 
With one iteration by Eq. (3.33), the initial guess f is updated to the value f , and the error 
function e is decreased. 

To find the optimum value of t, we solve the equation: 


de(T)_ 0 

dt " a 


(3.34) 


For example, if e(f) = (g - Kf) 2 + ylifj 2 , the optimum value of t becomes. 




a/; 


t = 




• v j 




(3.35) 


The physical significance of the yZfj 2 constraint will be discussed in the next chapter. 


3.2.4 Landweber iterations 

Landweber [1950] has proposed another iteration formula for solving Fredholm integral 
equations of the first kind: 

f N = f N ~ l +K{g-Kf). (3-36) 

Landweber has shown that if the kernel functions satisfy the constraint: \ \ K 2 (x,y,r') dx dy 
dV' < 2, then the residual e = (g - Kf) 2 must converge to zero and that the solution can be 

written as: 

f N = {KKY X [/-(/- KK) N+l ]Kg (3.37) 
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(see Chapter 7 of Twomey [1977] for a derivation of this result). In this expression, we have 
assumed that the initial guess f° = 0. Note that the constraint imposed on the kernels poses no 
additional problem, since we can always divide the system g = Kf by an appropriate scaling 
factor that is incorporated into the kernel elements. 

If N — » oo, Eq. (3.37) approaches an unstable result, f = (ATK) 'Kg (see section 4.1 for 

more on the unstable nature of this solution in the context of error magnification). Twomey 
[1977] has shown that if the number of iterations is kept small, however, some stability in the 
solution is retained. In fact, if N is not too large, Eq. (3.37) will filter small eigenvalues in 
much the same way as the constrained linear inversion technique with H = I (see section 3.2.1). 
Thus, a properly truncated Landweber iterative technique is a valid way of solving our 
problem. 

At this point, we can ask what is the relationship between the Landweber iteration 
method and the method of steepest descent? The connection can be made clear by simply 
letting t = 1/2 and all y’s = 0 in the method of steepest descent [see Eq. (3.33)]; the resulting 
iteration formula is then completely equivalent to the Landweber method [Eq. (3.36)]. [Note 

also that when e = (g - Kf) 2 , the gradient of e becomes: Ve = -2A*(g - Kf)]. 

Hence, the Landweber iterative method is a special case of the method of steepest 

descent. The benefit of the Landweber method, however, is that it converges to the absolute 

minimum of (g - Kf) 2 when the kernel magnitudes are properly scaled. This is a powerful 

result, especially in view of the fact that the Landweber iteration is basically equivalent to a 

gradient search procedure. By contrast, convergence to the absolute minimum of the error 

function given in Eq. (3.32) using the iteration formula given in Eq. (3.33), is not assured. 
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CHAPTER 4 


EXTERNAL CONSTRAINTS 

The primary concern of this chapter is to describe possible constraint functions, C,(f), 
C 2 (f), etc., that can be used in the method of steepest descent (see section 3.2.3). Since 

constraints tend to bias the solution toward specific types of distributions, it is important that 
the constraints be based on sound physical principles. The effects of some of the constraints 
presented here will be tested later on, in Chapter 5. We will also show that it is necessary to 
overcome biases in our charge solutions that are produced from variations in the kernel 
functions (see Figures 3.3 and 3 4). A kernel scaling procedure that is designed to remove such 
"kernel biases" will be discussed in detail. Later, in Chapter 5, it will be found that this 
procedure is needed to obtain acceptable solutions when the method of steepest descent, or 
Landweber algorithms are used. 

Before discussing these main points, however, we first examine in a more general 
context, why external constraints are necessary in solving our problem. In addition, since we 
will present Q- and P-model results in Chapter 6, we have also included a discussion of the 
constraining process involved in the Marquardt algorithm. 

4.1 The need for external constraints 

In Chapter 3, we introduced a linear procedure for inverting field changes to find an 
unknown volume source distribution. Our unknown distribution, f, is a (n x 1) column vector 
of source charges that produce m field change values on the ground. We store these field 
change values in a (m x 1) column vector g, and a (m x n) kernel matrix K relates f to g, i.e., g 
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= Kf. Of course, in any real problem there are errors in the measurements, so we may more 
accurately describe our problem by writing g + e = Kf, where e is a column vector of 
measurement errors. With this system of equations, there are three basic reasons for applying 
external constraints: (1) to reduce possible magnifications of error in AE, (2) to reduce solution 
ambiguities (i.e., to avoid nonphysical solutions), and (3) to remove kernel biases. 

The idea of error magnification has already been discussed briefly in section 2.3. In 
general (i.e., for arbitrary values of m and n), we see where error magnification arises by 

looking at a straight-forward "solution" to our linear system: f = (KK) 4 K(g +e) [Twomey, 
1977]. If KK has small eigenvalues (i.e., if KK is ill-conditioned), then the solution error 
given by (KK) 'e = adj(KK)Ke/det(KK) = adj(KK)Ke/TI^ becomes very large. Here, adj( ) 
and det( ) are the standard adjoint and determinant operators, respectively, and the X/s are the 
eigenvalues of KK. 

The only way to prevent such an error magnification is to effectively increase the small 
eigenvalues of KK. One can see from Eq. (3.31), that this is achieved by adding the constraint 

matrix yH to KK. This adds to each eigenvalue of KK a reasonably large number so that the 

new set of eigenvalues corresponding to the matrix (KK + yH) are all sufficiently large. The 
filtering of small eigenvalues is particularly obvious for the case H = I. Adding the 
eigen-equation for KK to the eigen-equation for I gives: 
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KKx = kx 
ylx = yx 


(4.1) 


(KAT + yl)x = (k+y) x 


The inversion of (KK + yl) involves new eigenvalues (k + y) that are all larger than k in 

magnitude and do not produce a large error magnification. 

The charge ambiguities that are inherent in this problem have been briefly discussed in 
Chapter 1 and Chapter 3. Again, one can uniquely determine the fields from a given set of 
charges, but to find the charges from a finite number of field observations is fundamentally 
ambiguous. For example, any spherically symmetric charge distribution with an outer radius R 
and total charge Q, will produce the same field at r > R , i.e., E = (Q/r 2 )r . If a cloud-to-ground 
lightning changes the cloud charge distribution in a spherically symmetric way, we will have 
difficulty in determining the true radial dimension R and charge distribution deposited. 

To remedy the situation, we can impose certain restrictions on our solution that are 
physically reasonable and are based on our knowledge of lightning phenomena. For instance, 
we might require that the charge density p(r) inside the thundercloud does not exceed some 
maximum limit; or we might restrict the radial dimensions of the charge distribution to be 
within, say, 4-5 km. 

The fact that more than one charge configuration can produce exactly the same field 
pattern on the ground can be demonstrated in a formal manner. If we split f into two parts, one 
part that is orthogonal to the row vectors of the kernel matrix K (defined by f 0 , where Kf c = 0) 

and one part that lies in the space spanned by the row vectors of K (i.e., a nonorthogonal part 
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f n ) we then have: g = Kf = K(f 0 + f n ) = Kf n (note that we are neglecting any measurement errors 

for this discussion). It can be shown (see Appendix) that the portion of the solution that lies in 
the row vector space of K, given by f n , is unique, i.e., this portion is totally constrained by the 

measurements. Conversely, the measurements are "blind" to variations in f Q (i.e., f Q is 
completely unconstrained by the measurements since Kf 0 = 0). The infinity of possible solu- 
tions is then a result of the infinity of possible choices for f 0 when m < n. The only way to 
constrain f 0 is to add external constraints (i.e., constraints other than the measurements). In 

general, these external constraints can be used to remove physically meaningless portions of 
both the orthogonal and nonorthogonal parts of the solution. 

Note also that there will be solution ambiguities due to measurement errors, since 
infinitely many solutions can be found to satisfy: g { - £j < EK^fj < g i + Note that this type of 

ambiguity is a common property of all linear inversion problems and is independent of the 
ambiguity described in the previous paragraph; the reader is urged not to confuse the two. The 
inequality above could also be written, with no loss in generality, by replacing fj with the 

nonorthogonal portion of the solution (fj) n ; this tells us that there is indeed a need to constrain 
f n as well as f c . 

Finally, we can add constraints to remove certain biases that are inherent in our kernel 
functions. If we examine Figures (3.3) and (3.4), we see that the method of steepest descent 
and Landweber algorithm will always tend to place more charge in regions where the 
magnitude of the kernel functions are large (see section 4.3.1). For instance, the 1/D 3 
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dependence in the kernels (D is the horizontal distance from a field mill site) tends to produce 
solutions that have greater amounts of charge over the network (i.e., at small D). There will 
also be a tendency to place more charge at lower altitudes when the source is over the network 
and more charge at higher altitudes when off the network, i.e., the kernel functions have 
maxima at altitudes z max = D/V2. 

4.2 The Marquardt algorithm as a constrained linear inversion result 
In Chapter 3 we have seen that the Marquardt algorithm is basically a combination of a 
gradient search and an expansion method for finding the minimum of a C 2 hypersurface. In 
order to illustrate the constraints that are inherent in this algorithm, we will show that it is, in 
fact, a form of a constrained linear inversion discussed in section 3.2.1 and in Chapter 6 of 
Twomey [1977]. The reader is also referred to Hoerl and Kennard [1970] and Marquardt 
[1970] for more detailed discussions of this point. Note that in these earlier publications, the 
term "ridge regression" is used instead of "constrained linear inversion." Once the constraining 
process is clearly understood, the advantages of the gradient-constrained linear inversion 
method over the Marquardt algorithm will become evident. 

We begin by defining a sensitivity matrix S of the form: 


5 s ^M' = 


“3M j 

3M \~ 

9ai 

da p 

dM m 

■ ■ I 

L"^r ••• 

da p _| 


(4.2) 
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where M'j a M/C; [note that a^-.ap, O; and M ; have all been defined in section 3.1.1 in the 
discussion of the Marquardt algorithm]. The change in each due to small changes in the 
parameters is: 

AM' = S Aa . (4.3) 

Obviously, we would like our model field changes to approach the lightning field change 
measurements AE. To do this, we choose AM' = AE' - M'. Since Eq. (4.3) is exactly of the 
form g = Kf, we can immediately write down an iteration formula to determine the parameter 
increments in a way that is analogous to the constrained linear inversion formula given in Eq. 
(3.31): 


Aa = (SS+yH) 1 S AM' . (4.4) 

But Eq. (4.4) is just the iteration formula for the Marquardt algorithm [Eq. (3.5)] when we let y 
= X and H = I (note here that SAM' = S(AE' - M') = V, and SS = A'). 

We conclude that the Marquardt algorithm is basically an iterative procedure for finding 
the minimum of an error function, e, of the form: 

e = (AM ' - S Aa) 2 + y(Aa) 2 , (4.5) 
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where AM' = AE' - M'. Hence, only the model increments Aa are being (externally) 
constrained in the Marquardt algorithm (i.e., the term yAa 2 forces the parameter increments to 
be small when yis large). Without directly constraining the model parameters, anomolous 
results are possible; this important point should be kept in mind when the Q- and P-results of 
Chapter 6 are reviewed. Note that in the gradient constrained linear inversion algorithm 
discussed in section 3.2.1, the model parameters are, in fact, directly constrained. In contrast, 
most investigators to date have only placed "weak" constraints on the range of model 
parameter values (e.g., charge altitudes were constrained to be positive for Q-model analyses 
[Jacobson and Krider, 1976]). 

Furthermore, if we consider a lightning source that is far from the field mill network, 
iterations by Eq. (3.5) will tend to make the (X,Y) model parameters large in magnitude (i.e., 
the model charges will become horizontally displaced from the network in the direction of the 

lightning source). This results in the elements of A' (= SS) getting smaller and smaller so that 

A' becomes increasingly ill-conditioned with each iteration of Eq. (3.5). We have previously 
shown that the addition of the term yl will remove this ill-conditioned nature. Unfortunately, 
the Marquardt algorithm decreases y while approaching the minimum in the C r 2 surface (i.e., 

upon locating the source with the model charges). Since (A'+yl) will also be ill-conditioned 
when yis small, and since illconditioned matrices are slow to invert, the net result is that there 
will be slow convergence toward sources lying well off the network. The large elements of 
(A’+yl) -1 will also amplify the errors in AE^ and possibly produce large parameter errors in the 

absence of sufficient external constraints. 
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4.3 Constraints used in the method of steepest descent 
We will now examine the constraints that we use in the method of steepest descent. We 
will also comment on the possible choices of the H matrix in the gradient constrained linear 
inversion method. 

4.3.1 Scaling constraint 

In section 4. 1 we mentioned that the kernel functions can introduce biases into the 
solution. This can be seen more clearly by writing out one component of the gradient term in 
Eq. (3.33): 

3/7 = ” 2 ? (&• ~ ? K ‘ifj) K * + Yi ^ + - • (4.6) 

Since f k is updated in proportion to 3e/9f k [see Eq. (3.33)] which is, in turn, propotional to K^, 

variations in will give similar variations in f k . One way to reduce such biases is to scale the 

elements of the kernel matrix so that the new set of kernel elements no longer have large 
variations. 

In order to scale our problem, we can pre- and post-multiply K by square diagonal 
matrices A and B. Now, Eq. (3.26) becomes: 
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(4.7) 


g = Kf 


g = KBf' 
Ag = AKBf' 

g' = K r 


where g = Ag, K' = AKB, and f = Bf. Equation (4.7) is now a new set of linear equations to 
be solved. 

The scaled kernel matrix elements can now be written as: = A^B^, but we still 

must find the m scalars A i} m a t , i = 1, -.m and the n scalars B^ = Pj, j = l.". n that optimally 
reduce the variability in the elements of K\ Initially, we selected oc’s and p’s that minimized a 
variance function of the form: 


* J 


(4.8) 


The critical points of this function are given by: 


a i =(Lp y ^)/(Ip y AT ( /) 2 

j J 

&=(LoL i K i j)/(La i K iJ ) 2 . 

i i 

The optimum scalars are then found by the following iterative technique: 

(1) Assume an initial guess of all a i5 

(2) Use Eq. (4.9) to find the optimum p j5 

(3) - tVF; {V = tyd/dad*-- 

(4) Return to step (2) until F is sufficiently small. 


(4.9) 
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With this technique we were able to reduce F by several orders of magnitude to a final 
value of 5.623E+04 for m = 25 and n = 4851. Unfortunately, the new set of scaled kernels still 
possessed undesirable biases. 

A better kernel scaling procedure is to set all cq = 1, and iterate each pj using P' - tV'F 

A A 

= Pi(d/c)Pi) + •■•+ P n (5/3p n )}, where the kernel smoothing function is given by: 


v (P) = S 2 k, h f + ® jfii-ft., at*. f 

+ (Py tfy-P;-*, Kij-n f + (PyATv-P/H, ) 2 ] (4.10) 

+ (py^-p;^^) 2 + (hKij-h^K^) 2 ] . 

Here we assume that there is a rectangular grid system with ti grid points in both the x and y 

directions (the number of grid points in the vertical is not important in the present discussion. 
In addition, the summations in Eq. (4.10) are performed only over the inner points of the grid 
system. The set of Pj’s that minimize produces a new set of kernels that are smooth, i.e., 

Pj^ij — Pj-iKjj.i (or equivalently K'^ = K^., $ (PjK^ - Pj.,!^.,) 2 is small; similar statements can 
be made for the other terms in X F. 

4.3.2 Maximum charge density constraint 

We have previously mentioned that one possible physical constraint is to set an upper 
limit on the charge density or the charge value at each grid point. Krehbiel [1986] has recently 
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reviewed the electrical structure of thunderclouds and reports charge densities that are between 
1 and 10 C/km 3 ; and larger charge densities might exist in certain localized regions of the cloud 
[Winn, Schwede, and Moore, 1974]. We can assume that comparable (order-of-magnitude) 
values of charge density might be the upper limits of those involved in lightning. 

To implement this constraint, we define a constraint function Cj (see section 3.2.3) as. 

Yi Ci ff) = Yi £/A (4<11) 

Note that the addition of each constraint term to the residual (g - Kf) 2 in Eq. (3.32) has the 
effect of changing the shape of our hypersurface e(f), and, in general, will change Ve(f). 

The effect of this constraint is to minimize the amount of charge it places at each grid 
point, and/or to spread out existing charge more evenly among the grid points. In the extreme 
case, Yi = °°> n0 charge will be allowed on the grid system. The details of how strongly a 

particular constraint should be weighted (i.e., what value of Y should be used) are discussed by 
Hoerl and Kennard [1970] and Twomey [1977], The values of Y we have used are given in 

Chapter 5. 

The constraint (4.1 1) can be easily included in the gradient constrained linear inversion 
method by setting the constraint matrix H = I (the identity matrix). (This is the same choice of 
H that was used in the Marquardt algorithm as discussed in section 4.2.) 
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4.3.3 Smoothing constraint 


A three-dimensional smoothing constraint that has a form similar to Eq. (4.10) can be 
written as: 

TiCi if) = y 2 ? [(fj—fj-i) 2 + (fj-fj+if 

J 

(JrSnf + (fj-fm'i 1 (4.12) 

This constraint prohibits large variations in the charge values between adjacent grid points. It is 
a three-dimensional analog of the first differencing constraint commonly used to smooth 
temperature profiles in satellite inversion problems [Twomey, 1977]. 

4.3.4 Focusing constraint 

This constraint, still being explored by the author at the time of this writing, is closely 
allied to the form of C,(f). In this constraint, an attempt is made to minimize the number of 

grid points that can have a significant amount of charge. We attempt to find solutions with as 
simple a charge structure as possible without becoming oversimplified. To implement this, we 
minimize a function that counts the number of grid points that have a significant charge. The 
constraint is: 

Y3C 3 (;f) = Y3 Zc y ;=Y3 , (4.13) 
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where 5 is a small positive constant. The magnitude of 8 determines the sensitivity of the 
counting procedure. For small values of 8 , a small charge on the j lh grid point will be counted 
because Cj = 1. If 8 is chosen large, Cj = 0, and the grid point will not be counted. Thus, with 

small values of 8 , the constraint inhibits the placement of charge on all grid points so that a 
huge increase in the error function given by y 3 L Cj = Y 3 £ 1 = Y3 11 does not resu lt- 


4.3.5 Conservation of charge constraint 

Conservation of charge is a fundamental law of physics that can be applied to our 
problem. Note that a cloud discharge, for example, can only move existing charge from one 
place to another or perhaps separate charges that previously existed in a neutral region. 

Overall, the discharge cannot create or destroy net charge. Hence, if we add up all changes in 
charge that were involved in a cloud discharge, the sum must be zero. Note here that we do not 
need to worry about charge transfers associated with convective or precipitation charging 
mechanisms that occur during the discharge process, since these effects have already been 
accounted for by using the time-varying Field correction described in Chapter 2. With this in 
mind, a conservation of charge constraint can be written as: 

Y<c< (/) = y< (Lfi) 2 • ( 4 - 14 ) 

j 

Ideally, if the cloud discharge was completely contained within our grid system and we 
had infinite grid resolution, the value of Y 4 should, in theory, be set to infinity in order to avoid 
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disobeying the law of charge conservation. However, as is the case with most inversion 
problems of this kind, it is generally not preferred to force the solution to a specific result using 
large weighting factors [Twomey, 1977]. 

A conservation of charge constraint can only be applied to cloud discharges because 
ground discharges remove net charge from the grid system in the U.H.S.. In the future, it may 
be possible to use this constraint to discriminate between ground and cloud discharges. For 
example, if we start with y4 reasonably large and find a small residual (g - Kf) 2 , then the source 
was probably a cloud discharge. 

The conservation of charge constraint can be implemented by setting H = 1 (the "one" 
matrix having as every element unity) into the gradient constrained linear inversion method. 
However, the inability of the (singular) one-matrix to remove small eigenvalues may lead to 
problems of error magnification. 

4.3.6 Other constraints 


There are many other constraints that could be used with our methods depending on the 
information that is available. For example, one could use weather radar data to help position 
the solution grid, and one could make high reflectivity regions within the grid more likely 
locations for lightning charge. Similarly, at the microphysical level, one could use 
measurements of cloud temperature, pressure, and moisture content to help position the 
solution grid, and to pre-bias regions within the grid that might be preferred locations for 
lightning initiation. 
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In summary, since there are an infinity of possible solutions to choose from, the 
addition of carefully selected constraints will allow us to find those solutions that are the most 
physically reasonable. If a constraint makes the residual (g - Kf ) 2 large, one should re-examine 
how realistic the constraint is, or the values of y used. 
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CHAPTER 5 


SOLUTION ERROR 

The problem of inferring the charge distributions that are deposited by a lightning, 
solely from ground-based measurements of AE has been reduced to solving a linear system of 
equations g = Kf. We have developed methods for solving this system, such as the method of 
steepest descent, and have found that a special case of this method, a Landweber iteration, has 
a number of advantages. Since there are inherent ambiguities in our solutions that arise from: 
(1) the non-uniqueness of charge and (2) errors in the values of AE, we have stressed the need 
for adding external constraints and have given some examples in Chapter 4. 

In this chapter, we will test the the method of steepest descent and Landweber 
algorithms using computer generated AE data. In these simulations, we first compute the fields 
from known lightning sources, add realistic measurement errors, and then apply the algorithms 
to find a charge solution. By comparing the results with the known lightning source, we 
compute solution error. 

Before presenting these simulations, however, we will first describe what a solution 
error is, how such an error is produced, and how the error is related to the independence of our 
measurements. The independence analysis will be given in the form of an eigenanalysis as 
described by Twomey [1977], Chapter 8. When we apply these methods to lightning data 
obtained at KSC, we will see that the geometry of the field mill network, the geometry and 
resolution of our source grid, the error in the AE values, and the number of measurements of 
AE £ all play a role in determining the information content of the field mill network. 
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5.1 Information content 


As a general rule, we expect that if we have more information in our measurements, we 
will obtain more accurate solutions for the lightning charge distribution. One can view the 
measurements as the primary or "intemal" constraints and any added constraint functions as 
"external". If the measurements are nearly redundant, our solution will not be clearly defined; 
in this case, we would say our intemal constraint is "weak" or that there is "little information" 
in the measurements. For example, if all 25 of the KSC field mill sensors were placed very 
close to one another, each sensor would detect almost the same field change, and the entire 
network would be roughly equivalent to a single measurement. In this case we would have 
very little information about the lightning events, and there would be little chance of finding 
accurate source distributions. 

5.1.1 The eigenanalysis test for independence 

The eigenanalysis test that we will use to determine the number of independent pieces 
of information that are contained in m measurements has been discussed by Twomey [1974] 
and Twomey [1977, Chapter 8]. Starting with our linear system, g = Kf, we write the 
unknown charge distribution in terms of a linear combination of orthonormal functions V;(r') 

(written discretely as v^ so that we have: 

m 

f = X a, v, = d>a . (5.1) 

i=i 
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Here, <I> is a (nxm) matrix whose m column vectors are Vj, i = l,...,m, and <J>C> = I. 

Orthonormality can be assured if we let d> = KUA-i' 2 where U is a (mxm) orthonormal matrix 

having as columns the m eigenvectors of the covariance matrix C m KK, and where A is a 

(mxm) diagonal matrix having as diagonal elements the m eigenvalues of C. Since any real 
symmetric matrix (such as C) can be written in terms of its eigenvalues and eigenvectors, i.e., 

C = UAU, we may verify the orthonormality of as follows: 

<DO = (A- 1/2 UK) (KUA 1 ' 2 ) 

= A- 112 U(KK)UA 1/2 

= A 1/2 (UCU)A* 1/2 (5.: 

= A-*/ 2 AA-w 
= 1 . 

We can now write g and f as follows: 

g = Kf = Kd>a = KKUA ,/2 a 

= CUA' I/2 a = UAA 1/2 a = UA 1/2 a (5.3) 

=> a = A 1/2 Ug 

/. f = KUA 1 Ug. 

The last line in Eq. (5.3) shows that if the measurements, g, have an error, e, a solution error, s 
= KUA-*Ue will result. 

The square magnitudes of f, g, e, and s are now: 
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f>=Za! 

I 

g J =ZX,a? (5.4) 

i 

e 2 = r 2 g 2 

F 12 

5 2 =Z^, e' = i/e, 

i A»i 


where r is a fractional number that is arbitrarily chosen to estimate the overall measurement 
error (i.e., r ■ e 2 /g 2 ). In the discussion preceding Eq. (4.1), we saw that large solution errors 
are generated when the measurement errors are divided by small eigenvalues. This result is 
apparent in the last line of Eq. (5.4). The eigenvalues depend solely on the geometry of the 
field mill network and on the geometry of the solution grid that is used to describe the lightning 
charge distribution. It is interesting to note that, given a particular network and grid geometry, 
each error component e i is amplified by a distinctly different eigenvalue. Thus the overall 

solution error is also sensitive to how the measurement errors are distributed across the 
network. 

Typical square magnitudes of the four vectors in (30) can be found by letting ^ = a^, = 
1/Vm (i.e., by scaling the problem so that f 2 = 1). This gives: 
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p= 1 


g 2 =ZXi(l/m) = <Xi> 

l 


(5.5) 


e 2 =r 2 <\> 

s 2 = — X J- = e 2 <l/Xi> = r 2 <X;> <1/Xi> , 
tn i Ki 


where the brackets o denote averages. 

We will now write down several criteria that we can use to determine the independence of 
m measurements. For a scaled problem with <X> - X m Jm and <\[k > ~ l/(mX min ), we may 

write: 


FORM 1 


If s 2 « f 2 = 1 (i.e., if » e 2 /m) 

then the independence of m measurements is assured. 


(F.l) 


or, equivalently: 


FORM 2 


then the independence of m measurements is assured. 


(F.2) 
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In general, the [1 th measurement within a set of m measurements will consist of three 
parts: (1) a predictable portion that is based on a linear combination of the other (m-1) measure- 
ments, (2) an unpredictable part that provides new information about the unknown distribution f j5 

and (3) an error term that is proportional to a linear combination of the measurement errors. 
These contributions can be summarized as follows: 


= £ bigi 

+ X kjfj — 

£ bi Ci , 

(5.6) 


j 

i 


predictable 

information 

error 


term 

term 

term 



where k (or k(r')) is the error in predicting the (3 th kernel function from all other kernels, i.e., 
ATp(r')= ZbiKii r') + k(r'). 

If the error term is large and/or if the kernel functions are highly dependent such that kj = 

0, the error term will tend to dominate the information term. In this case, the (3 th measurement is 
essentially dependent on the other (m-1) measurements (plus some error) and provides us with no 
new information about f. From the standpoint of our inversions, the p 1 * 1 measurement is not 

worthwhile. With this though in mind, we can write down another test for information content: 


FORM 3 


If (Ik/p 2 » (Ib^) 2 

then the independence of m measurements is assured. 


(F.3) 
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The equivalence of FORM 3 with FORM 1 and FORM 2 is given in Twomey [1977, Chapter 8]. 
In the eigenanalyses that we will give below, we will use FORM 2 of the eigentest. 

It should be noted that any set of m measurements are dependent if, and only if, (Ekjfj - 

ZbiEj) = 0. If each measurement within a given set provides new information about fj, then we 

can regard each of these measurements as being independent. If, however, each measurement 
does not provide new information about fj (i.e., if (Ib^) 2 » (Lk^) 2 ), then we will say that these 

measurements are "dependent" even though (Ikjfj - Zb^) * 0. 

5.1.2 Eigenanalyses of the KSC field mill network 

Before we can do an eigenanalysis of the KSC field mill network (see Figure 2.1), we 
must first select a grid system for the source. This is because the eigenvalues are derived from 
the covariance matrix C, whose elements depend on both the geometry of the field mill network, 
and the geometry of the grid in V'. (Note, a general discussion of the grid system has been given 
in section 3.2 and the volume V' is illustrated in Figure 3.2.) 

Obviously, we would like to have a large, high resolution grid system so that we can 
describe every discharge that occurs near the network. We would also like to have a grid that 
produces large eigenvalues so that we have small solution errors [see the last line in Eq. (5.4)]. 

To find an optimum grid, and to understand more about how the information content varies with 
different grid types, we have analyzed several grids with various grid sizes, resolutions, and 
locations (see below). 
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To find the eigenvalues of C, we have used a standard Jacobi transformation technique 
that has been described in Chapter 4 of Twomey [1977]. For each grid, the number of Jacobi 
transformations was typically 1400-1500. All eigenvalues and eigenvectors were checked by 

determining how closely UAU approximated C; we required agreement to 6 or more significant 

figures. 

In order to have a space-filling grid, we will begin with a horizontal grid dimension of 40 
x 40 km 2 , and a vertical dimension of 20 km. With these dimensions, we will try grid resolutions 
of 1, 2, and 4 km. Table 5.1 shows that the largest eigenvalues are obtained when the grid 
resolution is high, i.e., 1 km. Hence, there is a trade-off between improved resolution (reduced 
solution error) and computing time. However, the payoff in reducing the solution error by 
improving the grid resolution from 2 to 1 km is negligible (i.e., the relative error magnification, 
s 2 /! 2 , is on the order of 10' 2 in each case). (Note that the values of s 2 given here, and below, are 
directly comparable since all solution errors are relative to an f 2 value of unity.) 

To estimate the number of independent measurements that we have in the KSC network, 
we can assume a 10% random error in the values of AE (see section 2.3). This corresponds to a 
value of r = 0.1 in FORM 2 of the eigentest, or a critical value of (r/m) 2 = 0.16E-04, for m = 25. 
Many authors use a more conservative test (see Twomey [1977, Chapter 8]) and let m = 1, even 
though there are still 25 measurements. In this case, the critical value becomes (r/m) = 0.1E-01. 
Now, if we require the scaled eigenvalues, X-* = X. i A max , to be one order of magnitude larger than 

this conservative test value, we can see in Table 5.1 that we have 25, 10, and 5 independent 
pieces of information for the grids with 1 , 2, and 4 km resolutions, respectively. 
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Table 5.2 shows the eigenvalues of grids that do not fill the entire space but that preserve 
a one kilometer spatial resolution. Note here that there is more information content when the 
grid is large. In effect, it is more difficult to extract independent pieces of information about our 
unknown source distribution when this source is confined to a small volume of space. 

We will now examine the information content of the data as a function of grid location. 
To do this, we have computed the eigenvalues for a small cubic grid (dimension = 10 km, 
resolution = 1 km) at different grid locations. The results are summarized in Table 5.3. Note in 
this table that when the grid is located far from the network, less information is obtained. This 
expresses the fact that higher order moments of the unknown charge distribution are attenuated 
with distance. Since distant lightning sources produce similar fields at each field mill site, few 
independent pieces of information about the source are obtained. If the grid is displaced 
infinitely far away, C becomes singular, and the solution errors go to infinity (a case of zero 
information content). 

Using the above results, we have selected the grid system in column 2 of Table 5.1, i.e., 
dimensions 40 x 40 x 20 km , and a resolution of 2 km. This grid has reasonably large 
eigenvalues and allows us to adequately resolve a large volume of space without invoking large 
computation times. All analyses that follow (including the linear analyses presented later in this 
chapter and in Chapter 6) will be based on this grid. 

Now that we have selected a grid, we can investigate the effects of network geometry on 
the information content of the measurements. Generally, we expect that there will be small 
eigenvalues when the network covers a small area and large eigenvalues when a larger area is 
covered, provided that the network is not substantially larger than the area of our grid system. 
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Table 5.4 shows the eigenvalues of five different networks geometries: (1) a small-area network, 
KSC networks operated in (2) 1978, (3) 1987, (4) 1991, and (5) a large-area network. All 
networks had 25 field mills, except for the 1987 and 1991 KSC networks that operated 34 and 31 
mills, respectively. The small network was a square grid covering 5.76 km 2 (i.e., the X and Y 
locations of the mills ranged from 17 km to 19.4 km in steps of 0.6 km) and the large network 
covered 1024 km 2 (X and Y ranged from 1 km to 33 km in steps of 8 km). The smallest 
eigenvalue, A. min , for the large and small networks was 0.1721 and 0.8080E-09, respectively. If 

we use the conservative value of (r/m) 2 described above, i.e, r = 0. 1 and m = 1, we find that we 
have 25 independent pieces of information with the large network, and only 1 piece of 
information with the small network. By comparison, there are about 10, 11, and 1 1 pieces of 
information derived in the 1978, 1987, and 1991 KSC networks, respectively. Note that there are 
9 additional measurements in the 1987 network over the 1978 network, yet only 1 to 2 more 
pieces of information are being gained about the lightning charge distribution (similar comments 
hold when the 1991 network is compared with the 1978 network). As discussed above in 
relation to Eq. (5.6), it is not always assured that each additional measurement will provide 
additional information. 

Note that the inversion error magnification is greater with the 1987/1991 networks than 
with the 1978 network (i.e., the relative error magnification, s 2 /! 2 = s 2 /l = s 2 , for the 1978, 1987, 
and 1991 networks are 2.322E-02, 7.099E-02, 2.483E-02, respectively). But, the 1987 and 1991 
networks have at least one more piece of information than the 1978 network, as described above. 
This situation is not a contradiction however, since the eigentests described/applied above are 
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only used to determine the number of independent pieces of information from the total set of 


measurements. It is up to the inverter to determine which, from the total set, are independent. If 
this is not done, the mills that produce small eigenvalues will give rise to large error magnifica- 
tions. In effect, you can get 1 1 pieces of information from the 1987/1991 networks if and only if 
you pick the 1 1 or 12 measurements (of the 31 or 34 present) that are truely independent. If you 
allow extra (dependent) measurements into the inversion process, the overall information will be 
decreased due to excessive error magnification. Of course, an alternative to deleting redundant 
measurements is to filter the small eigenvalues (see section 4.1). 

In summary, given the errors in AE, not much additional information has been obtained 
by adding/moving mills (from 1978 to present). Instead, our results show that it is more 
desirable, albeit less practical, to spread out the 25 mills of the 1978 network into an orthogonal 
arrangement. Mills placed close to one another clearly give rise to serious error magnification in 
the absence of small eigenvalue filtering. 

5.2 Simulated lightning sources over the network 

In order to test how well the method of steepest descent and Landweber algorithms 
retrieve a volume charge distributions that is directly over the network, we have computed the 
field changes that would be produced by a known source and then have analyzed these fields 
with and without simulated measurement errors. We start by examining the Landweber iterative 
method [Eq. (3.36)] and illustrate the effects of kernel scaling. We have also tested the effects of 
some of the external constraints discussed in Chapter 4 using the method of steepest descent. 
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The solution errors that are associated with (point) sources lying off the network will be given in 
the next section. 

For our first test, we inverted the AE’s from a 40 C point source at: X = 18 km, Y = 18 
km, Z = 8 km. The results of the Landweber iterations are given in Figure 5.1. Note that here 
and in all solution plots given below, there is a 2 km grid resolution and that X and Y range from 
-2 km to 38 km and Z ranges from 0 km to 20 km. Contour units are indicated at the top of each 
plot, and the letters "H" and "L" in the figures designate regions of positive and negative charge, 
respectively. For our purposes, we are particularly interested in the shape of the inferred 
solution, the location of the charge centroid, and the total charge. 

The volume distribution that is summarized in Figure 5.1 involved 4851 grid points, 
solutions such as this typically required about 19 minutes of computing time on a VAX 
computer. In most cases, it took less than 1000 Landweber iterations to insure that the rms 
relative residual |g - K/|/Vm was smaller than the rms relative measurement error lel/Vm. It is 
important to emphasis that, because there exists errors in g, we do not try to iterate until each 
component of the residual vector, (g - Kf), is equal to zero. Once we know that the vector Kf is 
within e of g, no further improvement in the solution can be guaranteed. In fact, as discussed in 
section 3.2.4, excessive Landweber iterations that force the residual to zero result in spurious 
solutions. In the case of the method of steepest descent, the Twomey-Chahine method, or most 
other iterative methods where absolute convergence is not assured, it is usually difficult to make 
each component of the residual vector smaller than each component of e. At least one or two 
components of the residual vector are always found to be too large. 
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Note in Figure 5.1a, that the contours of iso-charge define almost perfectly concentric 
circles that are centered on the true (X,Y) source location. The altitude cross-sections in Figure 
5.1b, however, show excess charge near the lower boundary. Since the Landweber iteration 
updates charge at each grid point by an amount that is proportional to the weight of the kernel at 
that grid point (see section 4.3.1), we believe this excess charge is caused by kernel biasing. 

To illustrate the potential benefits of kernel scaling, we have applied Landweber 
iterations to the same point source as in Figure 5.1, but with a new set of scaled kernels. The 
values of the original kernels ranged from 0 to 4443, and the scaled values ranged from 0 to 
0.02774. The results with the scaled kernels are shown in Figure 5.2. Note in Figure 5.2 that the 
vertical cross sections have better symmetry and that the maximum charge is closer to the source 
height of 8 km. (A summary of the errors in these solutions and other results to be discussed 
below are given in Table 5.5.) 

For our next test, we analyzed the same point source as in Figures 5.1 and 5.2, but with a 
random 10% measurement error and a 30 V/m digitizing error added to each AE value (note: all 
remaining simulations include similar errors). The solution with random errors is summarized in 
Figure 5.3. Because truncated Landweber iterations effectively remove small eigenvalues (see 
section 3.2.4), we have avoided large random oscillations that are frequently present in linear 
solutions. Note that the symmetry of the solution in Figure 5.3 has been preserved. The centroid 
error is only 1.61 km, and the altitude of maximum charge is close to 8 km. The total charge, 
however, has been overestimated by 13 Coulombs and this corresponds to a 1.6 km 
overestimation in the centroid altitude. 
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Figure 5.4 shows the results that were derived for a centrally located (slanted) discrete 
dipole source, a typical intracloud discharge. The plan view at Z = 10 km shows that the lower 
positive charge at (X = 16 km, Y = 18 km) has been found; the upper negative charge, however, 
(Figure 5.4b) is closer to 7 km than to the true source height of 8 km. It appears that the lower 
positive charge has more influence on the solution than the upper negative charge. The lower 
positive charge is well pronounced, but the upper negative charge is lower in magnitude and has 
been "pushed” eastward by the positive charge. Nevenheless, the linear method has detected a 
substantial charge transfer across the center of the network in the presence of measurement 
errors. 

Figures 5.5 and 5.6 show some extreme cases: (1) a horizontal air discharge, and (2) a 
vertical intracloud flash, respectively. Figure 5.5 shows that the inversion has done reasonably 
well in retrieving both the negative and positive charge centers in the horizontal discharge. Since 
a large spatial gradient of charge is difficult to retrieve with smooth kernels, our results show a 
horizontal dipole of somewhat greater charge separation, but with a somewhat smaller charge 
magnitude, so that the charge moment is conserved. Dipole ambiguities of this type have been 
discussed previously by Kreibiel [1981]. In Figure 5.6 (vertical dipole known source), the lower 
positive charge dominates the measurements and Landweber iterations produce a solution that 
looks like a low altitude positive point source. The upper negative charge has been completely 
ignored in the solution. 

Next, we have attempted to retrieve a complex air discharge that is horizontal and is 
centered over the network. The discharge begins by bringing a total of 40 C from the north to 
the south, but then splits into two separate branches, one that deposits 20 C to the southeast, and 
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one that deposits 20 C to the southwest. Figure 5.7 describes the exact known source and 
provides our retrieval. Apparently the network is unable to resolve fine structure in the lightning 
source over small distances when the charges are of like sign, because the two charge centers to 
the south have been merged together into one large volume charge. The overall pattern, 
however, is reasonably good. 

We will now examine the effects of various constraint functions of the form, C(f), that 
were discussed in Chapter 4. In particular, we will test the effects of limiting the grid point 
charge, the smoothing constraint, and the conservation of charge constraint. These are 
preliminary tests, and of course, more work could be done to determine the optimum constraints 
for various types of lightning sources. 

Figure 5.8 shows the effects of limiting the charge density (i.e., of applying the maximum 
charge density or max p constraint discussed in section 4.3.2), when the same field changes 
associated with the known point source in Figure 5.3 are analyzed. Recall from Chapter 4 that 
the max p constraint is the constraint that attempts to minimize the charge values at each grid 
point so that p max is constrained to be small. In Table 5.5, we can see that there are large 

overestimations of the total charge when point sources were inverted. However, when the 
constraint y£jfj is used with y = 1.5xl0 8 (large values of yare used because the kernels are 

inversely related to ej, the total charge in the solution was reduced to 39.9 C, a mere 0.1 C from 

the known charge of 40 C. Unfortunately, this constraint also tends to spread the total charge out 
across more grid points. For instance, if there were only 2 grid points, and a total charge of I C, 
placement of all this charge onto one of the grid points would give Zf. 2 = (l) 2 = 1. However, if 
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1/2 C was placed at each grid point, we would obtain a smaller value Xf 2 = (1/2) 2 + (1/2) 2 - 1/4 

+ 1/4 = 1/2. This spreading out of charge has distorted the solution, and has resulted in a larger 
net centroid error equal to 2.30 km. 

Figure 5.9 shows that, in the absence of kernel scaling, we can obtain reasonable 
solutions for point sources by forcing the solution to be smooth. Using y = 5x10 in the 
smoothing constraint (see section 4.3.3), the solution has a maximum at the correct charge 
altitude (8 km). These results are in fair agreement with scaled kernel results given in Figure 5.3. 

Finally, Figure 5.10 shows the effect of the conservation of charge constraint (section 
4.3.5). Since this constraint should only be used for cloud discharges, we have tested it on the 
horizontal dipole given previously in Figure 5.5. The results are given in Figure 5. 10 with y = 

10 4 . Note that Figure 5.10 has three improvements over the less constrained solution given in 
Figure 5.5: (1) the location of the positive charge center has been improved; (2) the magnitudes 
of the positive and negative charges centered at Z = 10 km are slightly more equal (i.e., in Figure 
5.5a we have ILI/H = 249/138 = 1.80, compared to 273/154 = 1.77 in Figure 5.10a); and (3) the 
total charge (Xfj) is closer to zero [i.e., (Ifj) Fig 5 5 = 6.45 C, while (£fj) Kg . 5 .io = - 0.14 C]. The 

reader should note in Table 5.5, however, that the overall centroid error has increased from 0.333 
km to 0.96 km. 


5.3 Simulated point sources off the network 

We will now examine the solution errors for (point) source charges that are horizontally 
displaced from the measuring network. Since most storms occur off the network, these results 


67 



will help to determine at what range a discharge can be accurately located. 

Figure 5.1 1 shows the solution for a point source well to the east of the network when no 
kernel scaling is performed. Because the kernels have a 1/D 3 dependence, our inversion tends to 
place more charge closer to the network than far away. The position errors in Table 5.5 show a 
large error in X (X relrieved - X tnie = -6.25 km). When we invert the same source with scaled 

kernels (Figure 5.12), the error is improved to -3.34 km). 

Figures 5.12 through 5.15 show qualitatively similar results but for point sources to the 
south, west, and north of the network, respectively. All of the retrievals tend to be closer to the 
network than the true source. The errors in horizontal position are usually improved if kernel 
scaling is used. 
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CHAPTER 6 


STORM ANALYSIS 

In this chapter, we present Q- and P-model results for eight different storms at the NASA 
Kennedy Space Center in July and August 1978. We will also show the results of our linear 
inversion for three flashes and will compare these solutions with the associated Q- or P-model 
results. These comparisons will illustrate the limitations of both inversion methods. We 
conclude this chapter with a section on future work. 

6.1 Flashing rate histograms 

Figures 6. 1-6.4 show the flashing rate histograms corresponding to all lightning activity 
that was detected by the field mill network during the storm periods of interest. Each lightning 
event was identified manually by plotting the digital field mill data onto a high-resolution video 
monitor as discussed in Chapter 2, sections 2.2.1 and 2.2.6. Previous flashing rate histograms 
(e.g., Jacobson and Krider, [1976]; Livingston and Krider [1978]; Maier and Krider [1986]; and 
Koshak and Krider [1989]) were based on an automatic, less accurate method of lightning 
identification (see section 1). In addition, our manual method also eliminated a threshold 
criterion in the definition of a lightning event (usually 600 V/m at two or more sites) that tended 

to produce an underestimation of the true flashing rate. 

Figures 6.1 and 6.4 show that the storms on July 5, 6, and 31 are all large with peak flashing 
rates that exceed 60 flashes over a 5-minute interval. The other storms were small in comparison. 
In the following, we will compare the Q- and P-results for these small and large storms. 
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6.2 Q- and P-model results 


Figures 6.5-6. 1 1 show Q- and P-model results for all large lightning events, i.e., those 
lightnings that produced a AE ^ 1 kV/m at two or more sites. This selection criterion is 
somewhat less restrictive than the 1 kV/m at three or more sites used by Maier and Krider 
[1986], Koshak and Knder [1989], and Knder [1989], Each figure shows four plots: (1) the 
altitude of the Q- and P-solutions as a function of time; (2) plan views of the X-Y locations of the 
Q- and P-solutions relative to the field mill network; (3) projections of the Q- and P-solutions on 
a vertical plane that is oriented east and west; and (4) projections of the Q- and P-solutions on a 
vertical plane that is oriented north and south. Only plot types (1) and (2) are given for the 
storms on July 14 and August 13, 1978, because these storms were considerably off the network 
and have few optimum solutions. 

In each plot, the magnitudes of the Q-solutions are shown as circles with a radius r = 
(AQ/Ejj) , where represents an assumed value for the dielectric breakdown of the air [Koshak 

and Krider, 1989]. In all height vs. time plots, E b was 0.3 MV/m and in all plan views E b was 1.0 

MV/m, i.e., smaller circles in the plan views provided a less cluttered map. 

The P-vectors in all plots show the direction that positive charge has been effectively 
transferred by the discharge. Since we are analyzing changes in the cloud charge distribution, a 
P-vector that points from a negative change in charge to a positive change in charge is pointing in 
the direction of a positive charge flow (or current). In all height vs. time plots, the magnitude of 
the three-dimensional P-vector is plotted and each vector is rotated clockwise from the vertical 
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by the zenith angle 0 = arccos (AP Z AP)]. In the plan views and altitude cross sections, the actual 

projection of the three-dimensional P-vector is shown. 

Finally, since Table 5.5 shows solution errors to increase for distant sources and since 
Koshak and Krider [1989] found anomolous Q- and P-parameters for distant lightning, we have 
removed Q- and P-results associated with distant lightning from all altitude vs. time, and vertical 
cross-section plots given in Figures 6.5-6.11. Plan view plots show all solutions that were found 
within the boundaries of the KSC maps given. This allows the reader to see what solutions have 
been removed from the other plot orientations. 

The region of optimum accuracy of our Q- and P-parameters has been estimated by 
performing several simulated AE-inversions. Obviously, there are an infinity of possible source 
geometries and locations that could be studied in the simulations. To make things simple yet 
informative, we have only studied AE’s produced from known point charge sources. We expect 
that the error results will be generally similar to analyses that involve more complicated source 
distributions. 

To determine meaningful statistics of solution error, 100 Q-model inversions were 
performed at each (X,Y) location of the known point source. Each simulation had different 
simulated errors, but was always between 0-10%. The errors were also different from field mill 
to field mill. The known point source was always placed at an altitude of 8 km and was always 
assigned a charge of 25 C. The values of X and Y ranged from 0 to 35 km in 5 km steps so that a 
total of 64 source locations were studied for each storm, i.e., for each network geometry 
encountered. Since 100 inversions were done at each (X,Y) source location, a total of 6400 
inversions per network were completed. For known sources located over the network (e.g.. 
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X = 15 km, Y = 15 km), average position errors were always less than 0.5 km, and average 
charge errors were typically only a few tenths of Coulombs. Standard deviations were usually 
less than these average values. For sources located off the network (e.g., X = 30 km, Y = 30 
km), position and charge errors were larger (typically 1-2 km and 5-9 C, respectively). 

From our simulation results, all Q- and P-solutions that occur within the regions: 


STORM (1978) 

X RANGE (km) 

Y RANGE (km) 

July 5 

9-28 

3-33 

July 6 

11-29 

6-33 

July 11 

9-29 

3-33 

July 14 

9-29 

3-33 

July 17 

9-28 

3-33 

July 19 

9-29 

3-33 

July 31 

9-29 

3-33 

August 13 

9-28 

3-33 


are assumed to have position errors that are less than 1 km, and, in the case of Q-solutions, 
charge errors that are less than about 2 C. Any Q- or P-result lying outside their respective 
region of optimum accuracy given above is not plotted in Figures 6.5-6. 1 1 (except in the plan 
views). Again, this estimated region of optimum accuracy is based only on point source 
inversions having Q = 25 C and Z = 8 km, i.e., on typical cloud-to-ground charge parameters 
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[Koshak and Krider, 1989]. If smaller point charge sources were used in the simulations (e.g., 

Q = 1 C or 5 C) larger position errors, and proportionately larger charge errors would obviously 

result, thereby changing the optimum region of accuracy. 

Table 6.1 summarizes the total number of flashes that were identified in each storm, the 
number of flashes that were large enough to be analyzed based on our AE threshold requirement, 
the number of solutions that were obtained with a C r 2 < 10, and the number of optimum solutions 

that were found within the region of optimum accuracy. Note that we have only been able to 
model a small fraction of all flashes that were identified, and this should be kept in mind when 
we make inferences about the electrical structure of storms at KSC. 

The statistics of our optimum model parameters are summarized in Table 6.2. Note in 
this table that the two small storms on July 17 and July 19 have average Q-altitudes that are 
about 1-2 km lower than in the large storms. However, the small storm on July 1 1 has 
Q-altitudes that are similar to those found in the large storms. More detailed comparisons will be 

given later in this document. 

6.2.1 Storm on July 5, 1978 

As seen from the plan view plot in Figure 6.5, most of the lightning activity in this active 
storm occurred near the south end of the field mill network. A total of 1209 flashes occurred in 2 
hours and 20 minutes. 

One feature in the height vs. time plot in Figure 6.5 is the fact that all high-altitude 
P-vectors point downward, the low-altitude P- vectors point upward, and the mid-altitude 
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P-vectors tend to point horizontally. The Q-solutions are located within a narrow band from 
about 6 to 8 km and are co-located with the horizontally pointing P-vectors. This general pattern 
was found previously by Koshak and Krider [1989] and is a characteristic of all storms that we 
have analyzed. 

The west-east vertical projections in Figure 6.5 show that the low-altitude P-vectors tend 
to cluster below the Q-region. The high-altitude P-vectors, on the other hand, have a greater 
horizontal extent and are much more numerous. In addition, the west-east vertical projection 
shows two separate clusters of events that are separated by about 3 km. The cluster to the west 
has no low-altitude, upward P-vectors and was not as active as the main cluster to the east. 

6.2.2 Storm on July 6, 1978 

Figure 6.6 shows that this large storm was located northwest of the Field mill network. 
Koshak and Knder [1989] have previously analyzed a 20-minute portion of this storm that was 
centered on the time of peak activity, i.e., between about 1935 to 1955 GMT. Here, we have 
analyzed the entire life -cycle of the thunderstorm. Note that the pattern of Q- and P-solutions 
that we have obtained is very similar to the July 5 storm. 

The storm on July 6 began in the northwest and then moved to the southeast. In the 
height vs. time plot we show two positive Q-solutions (the hatch-shaded circles) that are sys- 
tematically higher than most negative Q-solutions. These postive events also occurred on the 
southeastern fringe of the storm, in the direction of storm movement. The largest postive 
solution occurred near the end of the storm at about 202230 GMT. All of these features are 
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reasonably consistent with what we might expect of a positive cloud-to-ground lightning event 
[Rust, MacGorman, and Arnold, 1981; Brook, Henderson, and Pyle, 1989]. 

The altitude cross sections on July 6 again show that the low-altitude P-vectors form a 
compact cluster under the main Q-solutions, and the high altitude P-vectors are more spatially 
extended. The abrupt termination of solutions to the west marks the westward extent of the 
region of optimum accuracy. 

6.2.3 Storm on July 11, 1978 

The height vs. time plots for this small storm (Figure 6.7) shows that the average 
magnitudes of P-vectors are smaller than those in the larger storms. Note that the scale has been 
change from the 600 C km in Figures 6.5 and 6.6 to 100 C km in Figure 6.7. It is also interesting 
to note in Figure 6.7 that there is only one low-altitude upward pointing P-vector and there is not 
a large separation between the high-altitude P-vectors and the Q-solutions as in the larger storms. 

6.2.4 Storm on July 17, 1978 

The results for this small storm, given in Figure 6.8, have many of the same 
charactersitics as the storm on July 11. Here, there are no upward pointing P-vectors at low 
altitude, and there is not a large separation between the high-altitude P-vectors and the 
Q-solutions. It is interesting to note, however, that there were two positive Q-solutions. The 
positive altitudes and magnitudes do not appear to be appreciably different than the normal 


75 


negative Q-solutions. The negative Q-solution occurring just before 1812 GMT has a very low 
altitude that may be the result of inadequate area filtering. 

6.2.5 Storm on July 19, 1978 

Note from Figure 6.9 that this small storm occurred in almost the same region as the large 
storm on July 5; hence, we can obtain a good comparison between the characteristics of large and 
small storms from these cases. Here we have far fewer low-altitude P-vectors pointing upward 
than the July 5 storm and very few optimum Q-solutions. One thing that is clear in Figure 6.9 is 
that the average magnitudes of the Q- and P-solutions are all smaller than in the July 5 storm. 

6.2.6 Storm on July 31, 1978 

A portion of this large storm has been analyzed previously by Koshak and Krider [1989], 
and our new results for the entire storm period are given in Figure 6.10. The distribution in 
Figure 6.4 shows that this storm had a bimodal flashing rate, and we have been able to extract 
accurate solutions (i.e., solutions with C r 2 < 10) only for the first portion of the histogram. A 

detailed analysis of the digital field mill records has shown that many flashes are overlapping in 
time during the second phase. This has made it very difficult to obtain accurate AE’s for 
individual discharges during this period, and we have obtained only a few acceptable solutions. 

In the height vs. time plot and the vertical cross sections in Figure 6.10, the Q-altitudes 
are much more variable than in Figures 6.5 and 6.6. In fact, many of the high- and low-altitude 
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P-vectors overlap the Q-solutions. Overall the Q- altitudes vary from just above 4 km to just 
below 12 km. 

6.2.7 Storms on July 14 and August 13, 1978 

These storms (shown in Figure 6.1 1) are included here only to illustrate the difficulties 
we have with distant storms. Most of the solutions having C r 2 < 10 also had X-Y locations lying 

outside the region of optimum accuracy. When the region of optimum accuracy was relaxed 
(i.e., the acceptable area was made larger), the additional solutions obtained were anomolous. 
Indeed, even with the current area filter used in the August 1 3 storm, many of the Q- and P- 
altitudes are anomalous; tightening up the area filter constraint would eventually eliminate all 
solutions from the plot. 


6.3 Linear method results and comparisons 

Now that we have seen the Q- and P-results for large and small storms, we will analyze 
three typical lightning events using the linear method. Two events had successful Q- or P- 
descriptions, and one event could not be described by Q- or P-models. 

Our first flash is taken from the large storm on July 6, 1978. This discharge occurred at 
201521.4 GMT and was successfully described by the P-model. The P-parameters were; X = 
14.6 km, Y = 23.1 km, Z = 10.3 km and AP X = 230 C km, AP y = -55.5 C km, AP Z = -123.15 C 

km. The linear method produced the distribution that is shown in Figure 6.12. Note that there is 
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a pronounced dipole-like charge distribution that is in the same general area of the network as 
our P-result. The centroid of the Ify-distribution is located at (19.2 km, 22.4 km, 8.3 km), and Xfj 

= 33.6 C. However, as we might expect for a complex flash, an extra region of positive charge 
was effectively deposited toward the southeast. Since asymmetries of this kind do not appear in 
our computer simulations of dipoles, it is possible that the discharge may have actually deposited 
positive charge in this region. We probably were able to obtain a satisfactory P-fit only because 
the positive charge to the southeast is not large. 

For our second example, we will analyze a discharge that occurred shortly after the flash 
discussed above (i.e., at about 202210.4 GMT). This event is of interest because according to 
our previous Q-results, it was a positive discharge to ground having parameters X = 17.3 km, 

Y = 17.1 km, Z = 10.7 km, AQ = -39.1 C (see Figure 6.6). The results of our linear inversion are 
shown in Figure 6.13, and are notably more complicated than a spherically symmetric charge 
distribution. A large amount of negative charge has been effectively deposited just east of the 
network, and this is consistent with positive charge being transferred to Earth in a positive 
ground flash. However, an extra amount of positive charge has been deposited to the southwest. 
The overall centroid of this distribution, (20.7 km, 22.7 km, 8.3 km), is still in fair agreement 
with the location of the Q-fit. And, the total charge change, Xf = - 41.6 C, is also very close to 

the Q-fit result given above. Since the negative charge (change) center dominates the solution in 
Figure 6.13, it is not difficult to understand why an acceptable Q-fit was found for this discharge. 

As our last example we will analyze one lightning event from the July 6 storm that we 
could not fit to either a Q- or a P-model. The results for this flash are shown in Figure 6. 14, and 
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it is clear why neither the Q- or P-model could describe this event. The plan view for Z = 10 km 
shows three prominent charge centers. We saw earlier, in Figure 6.12a, that a point-dipole model 
could describe three charge centers if one of the three charges is small, but this is not the case 
here. 


6.4 Future work 

In the future, we plan to make further analyses of natural lightning with the linear 
method, and will compare these results with Q- and P-model solutions. Thus far, the linear 
results are in fair agreement with the Q- and P-model results, and they also have explained why 
at least one lightning event failed to be described by a Q- or P-model. Further comparisons of 
this kind will help us to understand the limitations of both inversion procedures, and will eventu- 
ally lead to a better understanding of thundercloud electricity. 

In addition, note that all the work presented here on the new linear technique can be 
carried over to the study of a fundamentally identical inversion problem, where measurements 
of 9E/dt are made instead of AE, and where retrievals of 9p/9t are found instead of Ap. 
Specifically, one may write: 

E(x,y,t) = \ v .K(x,y,r') p(r',t)dV ' . (6.1) 

Taking a partial time derivative of this equation gives: 

^ (jc ,y,t) = \ v .K{x,y,r') 0 P ^ V) dV ' . (6.2) 
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This last equation can be solved in much the same way as the field change inversion problem 
discussed above. In effect, we are dealing with the same linear system g = Kf, except that the 
physical meaning of g ( dE/dx at the ground) and f (current, 3q/9t, in the cloud) have changed. 
Accordingly, the types of external constraints applied will, in general, be different from those 
discussed in Chapter 4. Since the kernel function in Eq. (6.2) is the same as before (i.e., the 
gradient of the Green function evaluated at the ground), all the same kernel scaling techniques 
discussed above can be applied here. Furthermore, the results of the information content eigen- 
analyses given above directly apply. Indeed, from a practical standpoint of obtaining reasonable 
inversion results, most attention is centered around successfully dealing with the kernel 
functions; the physical meanings of g and f are almost rudimentary. 

Using a ground-based network of dE/dl sensors that could also provide values of AE, or 
using field-inferred values of each quantity, our linear method could be used to help improve our 
understanding of the electrical nature of thunderstorms. Combinations of these inversion results 
with studies of the Maxwell current density at the ground (e.g., of the type described in Krider 
and Musser [1982]; Knder and Blakeslee [1985]; and Deaver and Krider [to be published in 
JGR]) should also be investigated. 

We can extend our above comments by noting that any temporal operator (e.g., d/dt, 

d 2 /d ^ , / ( )dt, etc.) can be applied to Eq. (6.1) and the same linear system g = Kf will result This 
system, in turn, can be solved by the methods presented in Chapter 3 above. Again, the types of 
external constraints will generally differ from one inversion problem to the next. 

Note that the inversion problem we are considering in Eq. (6.2) should not be confused 
with an interpolation algorithm discussed in Blakeslee [1984], for contouring current densities. 
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This algorithm, originally proposed by Dr. Elemer Magaziner (NOAA), cannot reconstruct 
realistic source distributions, but instead finds a single current pulse a fixed height above each 
field mill. The main purpose of the algorithm was to interpolate measurements of E (or 
measurements of AE, 3E/3t, etc.) for contouring purposes, and as such, is a far more restricted 
problem than the one we pose above. Furthermore, if more realistic current distributions are 
found by using the Landweber and/or steepest descent method described in Chapter 3, our 
methods would also serve as an improved interpolation scheme of E, 3E, 3E/3t, etc. As 
discussed in Blakeslee [1984], this interpolation is done by simply solving the forward problem 
after the source distribution is acquired. 

Finally, it is of great interest to supplement ground-based measurements of lightning AE’s 
with measurements of AE taken above ground with aircraft. Obviously, additional 
measurement(s) above ground will be highly independent of the ground-based network, so that 
the information content of the total ground and air measuring system will be much improved. 
This could ultimately lead to far more accurate lightning charge retrievals. 

The reader should note, however, that an aircraft measurement of lightning field change 
will in general be a vector quantity AE = (AE X , AE y , AE Z ). The ground-based measurements 

have, on the other hand, only a z-component value since, under our assumptions, AE X and AE y 

are always zero on the ground. Hence, before attempting an inversion, our basic Fredholm 
integral equation [Eq. (3.22)] must be further generalized. Taking the temporal change in the 
gradient of Eq. (3.17) leads to three linear Fredholm integral equations: 
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(6.3) 


gjr)= $ KJr,r')f(r')dV' 

V' 

g y (r)= $ K y (r,r')f(r')dV' 

V' 

g z (r)= jK 2 (r,r')f(r')dV', 

V' 

where the g’s are the changes in the x, y, and z-components of the potential gradient vector at 
location r, the K’s are the spatial gradients in the x, y, and z directions of the Green function 
given in Eq. (3.16), and f is the change in the charge density at r". These three equations can be 
added to obtain one Fredholm integral equation: 

g(r)=\K'(r,r')f(r')dV' . (6.4) 

V' 

Note that we have used a "g*" to indicate the sum of the field change components (g x + g y + g z ), 

and a "K*" for the sum of the kernel components (1^ + K y + K z ) to avoid confusion with Eq. 

(3.22). It is important to emphasize that, since the kernel K*(r,r') is different, and far more 
complicated than our old kernel K(x,y,r), the feasibility of inverting aircraft AE’s along with 
ground-based measurements is yet to be determined. In particular, the scaling of K*(r,r') is 
complicated and requires additional investigation. At present, it appears that it is most practical 
(both mathematically and experimentally) to invert only the last (z-component) equation in (6.3), 
i.e., to invert ground-based measurements together with any number of aircraft measurements of 
AE Z . Simulated inversions of this equation are currently being investigated by the author. 
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CHAPTER 7 


SUMMARY 

We have developed an interactive computer program that plots digital field mill data and 
can be used to identify and compute the changes in the cloud electric field due to lightning. We 
have then used this program to compute accurate lightning AE-values and flashing rates in eight 
Florida thunderstorms at the NASA Kennedy Space Center. The lightning field changes have 
then been analyzed using a nonlinear least-squares minimization procedure to derive the 
optimum parameters for point charge (Q) and point dipole (P) models. 

Three storms analyzed had relatively low flashing rates, while three other storms were 
considerably more active. We have also illustrated the difficulty in analyzing two distant storms. 
Because of the long time periods of storm activity studied (e.g., over several hours in most 
cases), the improved detection and analyses of lightning AE s, and reduced threshold biases, we 
have been able to make an improved comparison between large and small storms. Even with 
these differences however, the spatial and temporal development of the Q- and P-solutions were 
found to be similar to results given in Maier and Knder [1986], Koshak and Knder [1989], and 
Krider [1989]. 

We have also developed a new way to analyze lightning field changes. The new 
approach starts by describing the inversion problem in terms of a linear system of equations 
given by; g = Kf, where g is a column vector of field change measurements, f is a column vector 
that specifies the values of lightning charge on a grid system lying above the measuring network, 
and K is a kernel matrix relating the two. 
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We have reviewed some standard linear inversion methods, but because the kernel matrix 
K is very large, and since positive and negative charges must be retrievable, we have been forced 
to devise a new technique of inversion. The new (iterative) technique has been strongly 
motivated by the Twomey-Chahine tomographic inversion method and is essentially based on the 
method of steepest descent. Each charge value is updated via: f/ = fj - tfe/tJfj, where e is an 

error function of the form: e(f) = (g - Kf) 2 + Y,C,(f) +...+ Y N C N (f), and the C’s are constraint 

functions. A distinct benefit of this approach is that one can select any type of external 
constramt(s) desired, i.e., the constraint function are arbitrary. 

By letting t = 1/2 and removing the constraint terms C,(f),...,C N (f), we find that the 

iterative method produces a stable solution provided the number of iterations is limited (e.g., we 
have typically used 1000 iterations or less in our analyses presented above). We have shown that 
this special case of the steepest descent method is equivalent to an iteration technique introduced 
by Landweber [1950], 

Our linear method allows us to determine volume distributions of lightning charge on a 
grid of any desired resolution. To test the method, we have performed simulated inversions that 
use computer generated AE-values produced from known charge sources. In most cases, both the 
symmetry and (centroid) location of the retrievals were found to be in fair agreement with the 
known source. For sources located off the network, systematic errors in source location have 
been quantified and explained. We feel that the new linear method should be further explored to 
determine solution errors for source locations and geometries that differ from those presented 
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here. We also strongly recommend that additional external constraints be explored (see, for 

instance, section 4.3.4) to further reduce solution ambiguities. 

While we were developing the new linear method, we clarified some fundamental 
problems in field change inversions. We discussed and distinguished the difference between 
solution ambiguities that result from the nonuniqueness of charge and those that result from 
measurement error. And, we have clarified how it is possible to find lightning solutions when 
the number of unknowns is greater than the number of measurements made. We have pointed 
out that, as early as Workman, Holzer, and Pelsor [1942], the number of charge model 
parameters used to describe a lightning event have been (incorrectly) limited by the number of 
measurements. In our work, we have shown that it is possible to obtain valid results of up to 
n = 485 1 unknowns using just 25 measurements, provided external constraints are added. 

Moreover, even without applying external constraints to the model parameters, we have 
shown that it is still possible to solve for more unknowns than measurements. In the Marquardt 
algorithm, we have emphasized that the conditioned approximative curvature matrix A = A' + Y 1 
is invertible for fewer measurements than model parameters because the parameter incremen ts 

are constrained to be small. 
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APPENDIX 


UNIQUENESS OF THE SOLUTION LYING IN THE ROW SPACE OF K 

In Chapter 4, we have written the linear system of equations as: g = Kf = Kf 0 + Kf n , 

where f Q and f n are the orthogonal and nonorthogonal portions of the solution, respectively. 

Methods of inversion that restrict the solution to the space spanned by the row vectors of K are 
essentially forcing f c = 0. With this restriction, we will now show that the solution f n is unique. 

Note that for this discussion we will assume no measurement errors, i.e. the components of g are 
known exactly. 

Proof (by contradiction): 

Given: g = Kf 

Assume that there are solutions of the form f = la^ = K a where is the i** 1 row vector of K. 

In effect, we assume that f 0 = 0. From this set of solutions, pick two distinct solutions: 

r, = 

f 2 = l?a 2 . 
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=> Kf, = Kf 2 
=* K(fj - f 2 ) = 0 


=> KK(a, - a 2 ) = 0 

=> bKKb = 0, where bE(a,-a 2 ) 


=» (Kb) 2 = 0 
Kb = 0 . 


But since the row vectors of K are assumed to be linearly independent, Kb - 0 has only the 


trivial solution b = 0, so that f 2 = f 2 . Q.E.D. 
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Figure 2.1 Map of the KSC field mill network used in 1978. 


92 













Figure 2.2 A typical electric field mill site 
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Figure 2.3 An example of KSC digital field mill data. 
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Figure 2.4 An example of an F-change structure (mill #2, second flash). 
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Figure 2.5 An example of a more complicated F-change (mill #2, second flash) 
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Figure 3.1 Geometry associated with the a) Q-model and b) P-model. 
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Figure 3.2 Geometry used in the Fredholm integral formulation. 
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Figure 3.4 Kernel function (times 27T£.) vs. altitude, for various distances from 
f ieldmill . 




Figure 3.5 Schematic diagram of the method of steepest descent 
used to invert the linear system: g = Kf; note that Ai = -tVe, 
where t is an adjustable parameter. 
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Contours in 0.1 millicoulombs Contours in 0.1 millicoulombg 
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Figure 5.1 a Plan views of solution when no kernel scaling is per formed. 
Known source: ( 1 8Km, 1 8km , 8Km, 40C). 
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Figure 5.1 b Altitude cross-sections of solution for the source described, 
in previous figure. 





Figure 5.2a Plan views of solution when kernel scaling is performed. 
Known source: ( 1 8Km, 1 8km, 8Km, 40 C). 
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Figure 5.2 b Altitude cross-sections of solution for the source described 





Contours in millicoulombs Contours in 0.1 millicoulcmbs 
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Figure 5.3a Plan, views of solution when kernel scaling is performed, and 
measurement error is added. Known source: (1 8Km. 1 8km,8Km,40C). 
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Fiqure 5.3b Altitude cross-sections of solution for the source described 
in previozcs figure. 
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Figure 5.4 a Plan views of solution. Known source : ( 1 GKm, 1 8km, 8Km, 40C ), 
and ( 20 Km, 1 8 Km, 1 2Km,-40C). 
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Figure 5.5a Plan views of solution. Known source: ( 1 6Km, 1 8km, 8Km,40C), 
and ( 20 Km, 1 8 Km, 8 Km, —40C). 
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Figure 5.5 b Altitude cross— sections of solution for the source described 
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Figure 5. 6a Plan views of solution. Known source: ( 1 8Km, 1 8km, 8Km,40C), 
and ( 1 8 Km, 1 8 Km, 1 ZKm, — 40C). 
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Figure 5.6b Altitude cross— sections of solution for the source described 
in previous figure. 
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Figure 5. 7a Plan views of solution. Known source: (1 6 Km, 1 GKm, 1 0Krn,20C ), 
(20Km, 1 6Km. 1 0Km,20C), and ( 1 8Km,20Km. 1 0Km, —40C). 
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Fiqure 5. 7b Attitude cross— sectioizs of solution for the source described 
in previous figure. 
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Figure 5.8 a Plan views of solution. Known source: (1 8Km, 1 8km,8Km,40C). 
Maximum charge density constraint added. 
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previous figure. 





Contours in millicoulotnbs Contours in 0.1 millicoulombs 
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Figure 5.9a Plan views of solution. Known source: (1 8Km, 1 8 km, 8 Km, 40 C) 
Kernel scaling replaced by smoothing constraint. 
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Figure 5.9b Altitude cross-sections of solution; source describedin 
previous figure. 
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Figure 5. 1 0 a Plan views of solution. Known source: (16Km, 1 8km, 8Km, 40C), 
(20Km, 1 8Km,8Km, -40C). Conservation of charge constraint added. 
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Fiqure 5. 10b Altitude cross— sections of solution; source described in 
preirious figure. 





Contours in 0.1 millicoulombs Contours in 0.01 millicoulontos 
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No kernel scaling performed. 
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Figure 5.11b Altitude cross-sections of solution f or the source described 
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Figure 5. 12b Altitude cross-sections of solution for the source described 
in previous figure. 
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Kernel scaling performed. 
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Figure 5. 13b Altitude cross-sections o f solution for the source described. 
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Figure 5.1 4a Plan views of solution. Known source: ( 4Km, 1 8Km, 8Km,40C). 
Kernel scaling performed. 
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Figure 5.1 5a. Plan views of solution. Known source: ( 1 8Km, 32Km, 8Km, 40C ). 
Kernel scaling performed. 





Contours in millicoulombs Contours in millicoulombs 


I— I- 4 -I 
L_ J — L -I 
LL1J 

I I I I 


r r t i 

r r t i 

4— 4- 4 -4 
,LLi J 


r r t i 

hftl 
1-4-4 4 
LL1 J 
till 


— I — !— I- -I- 4 
LU 
_l_l_ L LI 
l I I I I 


I I I I 

r r t “i 

i- t- ■+ -t 
(-4-4 4 


\ I 

(r r r t 
I- r t t 
1-1-4 4 


*■"4 4 7 + 

- i- yj 

. l joi 
1 '0l 

TT-l 


4- K> 4 


4 I I A 

r t i/~] 

t-rxp 

I — 4 4*4 

llDj 

I I 1 I 


L i_ l j 

1 1 \ 1 
r r tv "i 

I- i- -tVt 

1-444 


Figure 5.1 5b Altitude cross-sections of solution for the source described 
in previous figure. 






Figure 6.1 Five minute flashing rate histograms for July 5 and July 6, 
1978. 
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6.4 Five minute flashing rate histograms for 
13. 1978. 
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Height vs. time and plan views. 
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Altitude cross-section views 
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Figure 6.6 Q— and P— results for storm on July 6, 1978: 
Height vs. time and plan views. 
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Altitude cross-section views. 
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Figure 6.7 Q- and P- results for storm on July 11, 1978: 
Height vs. time and plan views. 
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Altitude cross-section views. 
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Figure 6.8 Q- and P- results for storm on July 17, 1978: 
Height vs. time and plan views. 
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Altitude cross-section views. 
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Figure 6.9 Q— and P- results for storm on July 19, 1978: 
Height vs. time and plan views. 
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Altitude cross-section views. 
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Figure 6.10 Q- and P— results for storm on July 31, 1978: 
Height vs. time and plan views. 





147 


Altitude cross-section views. 
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Figure 6.11 Q- and P- results for two distant storms on 
July 14 and August 13, 1978: Height vs. time and plan views. 
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Figure 6.11 cont 
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Figure 6.12 a Plan views of solution for lightning flash, occurring at 
201 521 4GMT on July 6, 1978. Kernel scaling performed. 
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Figure 6. 12b Altitude cross-sections of solution for the source described 
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Figure 6. 13a Plan views of solution for lightning flash occurring at 
202210.4GMTonJuly6, 1978. Kernel scaling performed. 
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Figure 6.13b Altitude cross-sections of solution for the source described 
in previous figure. 





Contours in 0.1 millicoulombs Contours in 0.01 millicouloonbs 
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Figure 6. 14b Altitude cross— sections of solution for the source described 
in previous figure. 





TABLE 5.1 Eigenanalysis of a centrally located grid system 
for various grid resolutions. The center of the grid is 
always at (x,y,z) = (18 km, 18 km, 10 km) and has horizontal 
and vertical dimensions of 40 km and 20 km, respectively. 
Values in the table are the scaled eigenvalues A . s = k./k 
Expected solution errors s 2 are also shown for each grid. max 



Grid 

i 

Resolution (km) 

2 4 

01 

0. 1000E+01 

0.1000E+01 

0. 1000E+01 

02 

0.5Q2QE+00 

0.4196E+00 

0.3360E+00 

03 

0.4167E+00 

0.3177E+00 

0.2157E+00 

04 

0.3953E+00 

Q.2839E+00 

0.178QE+00 

05 

0.313QE+00 

0.2024E+00 

0. 1 107E+00 

06 

Q.2940E+00 

0. 1863E+00 

0.9028E-01 

07 

0.2633E+00 

0. 1463E+00 

0.5914E-01 

08 

0.2551E+00 

0. 1382E+00 

0.5227E-01 

09 

0.2388E+00 

0. 1227E+00 

0.4419E-01 

10 

0.2182E+00 

0.1049E+00 

0 . 3363E-0 1 

11 

0.2153E+00 

0.9591E-01 

0.2876E-01 

12 

0.21 14E+00 

0.9308E-01 

0.2380E-01 

13 

0. 1912E+00 

0.7782E-01 

0. 1900E-01 

14 

0. 1886E+00 

0.7434E-01 

0. 1339E-01 

15 

0. 1805E+00 

0.6560E-01 

0. 1307E-01 

16 

0. 1775E+00 

0.621 IE-01 

0. 1246E-01 

17 

0. 1718E+00 

0.5945E-01 

0. 1088E-01 

18 

0.1647E+00 

0.5396E-01 

0.8755E-02 

19 

0. 1570E+00 

0.4834E-01 

0.6600E— 02 

20 

0. 1488E+00 

0.4450E-01 

0.4447E-02 

21 

0.1476E+00 

0.4049E-01 

0.3157E-02 

22 

0. 1431E+00 

0.3314E-01 

0. 1396E-02 

23 

0. 1298E+00 

0.2949E-01 

0.8359E-O3 

24 

0. 1250E+00 

0.2544E-01 

0.5052E-03 

25 

0.1175E+00 

0.2005E-01 

0.1398E— 03 

S 2 

1 .274E-02 

2.322E-02 

4.508E-01 
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TABLE 5.2 Eigenanalysis of a centrally located grid system 

for various grid dimensions. Grid resolution is always 1 km, 

and the center of the grid is always at (x,y,z) = (18 km, s 18 

km, 10 km) . Values in the table are scaled eigenvalues A >i = 

A./X . Expected solution errors § 2 are also shown for each 
i • ..max 
grid. 



Cubical Grid 
5 10 

Dimension (km) 
15 

20 

01 

0. 1000E+01 

0. 1000E+01 

0. 1000E+01 

0. 1 000E+01 

02 

0.2077E-01 

0 .91 5 1 E — 01 

0.251 1E+00 

0 .4755E+00 

03 

0. 1827E-01 

0.8307E-01 

0.2205E+00 

0.4272E+00 

04 

0.6230E-02 

0.2368E-01 

0.1 141E+00 

0.3476E+00 

05 

0.8597E-03 

0. 1649E-01 

0. 1001E+00 

0.3162E+00 

06 

0.3419E-03 

0.8046E-02 

0.6768E-01 

0.2866E+00 

07 

0.3121E-03 

0.4469E-02 

0.5391E-01 

0.2831E+00 

08 

0 . 1 096E-03 

0.2913E-02 

0.4942E-01 

0.2562E+00 

09 

0.8020E-04 

0.2526E-02 

0.3265E-01 

0 .2379E+00 

10 

0.1800E-04 

0. 1 194E-02 

0.2402E-01 

0 .2131E+00 

11 

0. 1 130E-04 

0.9595E-03 

0. 1850E-01 

0.2083E+00 

12 

0.5051E-05 

0.5321E-03 

0. 1659E-01 

0. 2055E+00 

13 

0.3460E-05 

0.3358E-03 

0. 1452E-01 

0. 1968E+00 

14 

0. 1796E-05 

0. 1889E-03 

0.1194E-01 

0. 1898E+00 

15 

0.7817E-06 

0. 1496E-03 

0. 1040E-01 

0.1 766E+00 

16 

0.6453E-06 

0. 1409E-03 

0.8702E-02 

0. 1 723E+Q0 

17 

0.3875E-06 

0.9240E-04 

0.7818E-02 

0. 1661E+00 

18 

0 . 1 669E-06 

0 . 4667E-04 

0.6484E-02 

0. 1624E+00 

19 

0 . 7287E-07 

0.321 IE-04 

0 . 3769E-02 

0. 1495E+00 

20 

0.5334E-07 

0 . 1 779E-04 

0. 1864E-02 

0. 1438E+00 

21 

0.2123E-07 

0.141 8E-04 

0. 1034E-02 

0. 1318E+00 

22 

0.7835E-08 

0.5724E-05 

0.5089E-03 

0.2064E-01 

23 

0. 1590E-08 

0.2085E-05 

0.4213E— 03 

0. 1477E-01 

24 

0. 1014E-08 

0.4923E-06 

0.7132E-04 

0.3763E-02 

25 

0.4784E-09 

0.4269E-06 

0.4786E-04 

0 . 1447E-02 

S 2 

6.576E+04 

1 .038E+02 

1.353E+00 

1.083E-01 
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TABLE 5.3 Eigenanalysis of a small cubical grid system 
(dimension 10 km, resolution 1 km) for various distances from 
the measuring network. Center of grid is always at an 
altitude of 8 km. Values in table are scaled eigenvalues A. s 
= A. . For very illconditioned cases (i.e., when grid is 
well orr network) , the smallest eigenvalues could not be 
accurately computed and hence have been omitted below. 
Expected solution errors are given only for those cases where 
all eigenvalues have been found. 



(x,y) 

(18,18) 

Location of 
(25,25) 

Grid System 

(32,32) 

Center 

(39,39) 

(km) 

(46,46) 

01 

0. 1000E+01 

0. 1000E+01 

0.1000E+01 

0. 1000E+01 

0. 1Q00E+01 

02 

0 . 1 722E+00 

0. 1755E+00 

0.2592E-01 

0 . 3983E-02 

0. 1568E-02 

03 

0. 1579E+00 

0. 1 101E+00 

0. 1648E-01 

0.3402E-02 

0.9226E— 03 

04 

0.5838E-01 

0.4331E-01 

0.1041E-02 

0.51 26E-04 

0.5855E-05 

05 

0.4414E-01 

0.2228E-01 

0 . 7562E-03 

0 . 2580E-04 

0.2104E-05 

06 

0.3520E-01 

0. 1712E-01 

0.5174E-03 

0.8809E-05 

0.8335E-06 

07 

0. 1683E-01 

0. 1226E-01 

0.335 IE-04 

0.2689E— 06 

0. 1340E-07 

08 

0.1 531 E — 01 

0 . 8082E-02 

0.2154E-04 

0. 1872E— 06 

0.4766E-08 

09 

0. 1208E— 01 

0 . 3632E-02 

0.11 99E-04 

0.8320E-07 

0.3026E-08 

10 

0.7505E-02 

0. 1443E-02 

0 . 7062E-05 

0.3800E-07 

0.7670E-09 

11 

0.5937E-02 

0. 1155E-02 

0. 1165E-05 

0.2435E-08 

0.3181E-09 

12 

0.4458E-02 

0.485 IE-03 

0.4494E-06 

0.1141E-08 

0. 1529E-09 

13 

0.2549E-02 

0.3170E-03 

0. 1896E-06 

0.3151E-09 

0.1284E-09 

14 

0. 1383E-02 

0. 1274E-03 

0.6925E-07 

0.2504E-09 

0. 1057E-09 

15 

0.8865E-03 

0.451 IE-04 

0.5212E-07 

0 . 1 770E-09 

0.7397E-10 

16 

0.4913E-03 

0 . 2804E-04 

0. 1095E-07 

0. 1245E-09 

0.4279E-10 

17 

0.2685E-03 

0. 1552E-04 

0 . 7298E-08 

0.7937E— 10 

0.2903E-10 

18 

0. 1475E-03 

0 . 1 005E-04 

0 . 1 469E-08 

0.3725E-10 

0.1 361 E— 1 0 

19 

0. 1216E-03 

0.6217E-05 

0.4730E-09 

0.9219E— 11 

0. 1948E-1 1 

20 

0.7002E-04 

0.2158E-05 

0.3384E-09 

0.2854E-1 1 


21 

0 . 2824E-04 

0.6079E-06 

0. 1263E-09 



22 

0. 1492E-04 

0.3107E-06 

0 . 6439E— 1 0 



23 

0. 1004E-04 

0. 1649E-06 

0.2914E— 10 



24 

0. 1 102E-O5 

0. 1743E-07 

0.891 IE-11 



25 

0.6529E-06 

0 . 5484E-08 




S 2 

0.6583E+04 

0.5617E+06 





158 



TABLE 5.4 Eigenanalyses of five different field 2 
works: (1) a small square network of f rea . 5 *7. 6 . H;? 

1978 KSC network, (3) the 1987 KSC network, (4) the 1991 KSC 
network, and (5) a large square network of area 1024 km i . 
The U.H.S. grid has 2 km resolution with x and y ranging from 
-2 to 3 8 km, and z from 0 to 20 km. Expected solution errors 
s 2 are shown for each network. 


Small Area KSC KSC 

Network Network Network 

(1978) (1987) 


KSC Large Area 

Network Network 

(1991) 


01 

0.1000E+01 

0. 1000E+01 

02 

0.6410E-01 

0.4196E+00 

03 

0.6393E-01 

0.3177E+00 

04 

0.9410E-02 

0 . 2839E+00 

05 

0.7557E-02 

0.2024E+00 

06 

0.6190E-02 

0. 1863E+00 

07 

0.1 25 1 E — 02 

0. 1463E+00 

08 

0. 1086E-02 

0.1382E+00 

09 

0.2477E-03 

0. 1227E+00 

10 

0. 1864E-04 

0.1049E+00 

11 

0.1784E-04 

0.9591E-01 

12 

0.3274E-05 

0.9308E-01 

13 

0. 1376E-05 

0.7782E-01 

14 

0.9867E-06 

0.7434E-01 

15 

0.8476E-06 

0.6560E-01 

16 

0.3069E-06 

0.621 IE-01 

17 

0.2147E-06 

0.5945E-01 

18 

0. 1027E-06 

0.5396E-01 

19 

0.7777E-07 

0.4834E-01 

20 

0 . 2496E-07 

0.4450E-01 

21 

0. 1040E-07 

0.4049E-01 

22 

0.9157E-08 

0.3314E-01 

23 

0.1842E-08 

0.2949E-01 

24 

0.8621E-09 

0.2544E-01 

25 

0 . 8080E-09 

0.2005E-01 

26 



27 



28 



29 



30 



31 



32 



33 



34 




0.1000E+01 

0. 1000E+01 

0.5277E+00 

0.5202E+00 

0.3702E+00 

0.3587E+00 

0.2946E+00 

0.2844E+00 

0 . 2309E+00 

0.2065E+00 

0.2154E+00 

0.1914E+00 

0. 1661E+00 

0.1608E+00 

0.1491E+00 

0.1454E+00 

0.1341E+00 

0.1320E+00 

0 . 1 229E+00 

0.1155E+00 

0.1056E+00 

0.1122E+00 

0.9682E-01 

0.9389E-01 

0.8085E-01 

0.8822E-01 

0.7964E-01 

0.8275E-01 

0.7535E-01 

0.7984E-01 

0.7320E-01 

0.7717E-01 

0.6987E-01 

0.7129E-01 

0.6355E-01 

0.6755E-01 

0.6098E-01 

0.5681E-01 

0.5150E-01 

0.5503E-01 

0.4938E-01 

0.5153E-01 

0.4794E-01 

0.4444E-01 

0.4315E-01 

0.4275E-01 

0.3677E-01 

0.3759E-01 

0.3417E-01 

0.3475E-01 

0.3252E-01 

0.3189E-01 

0.2837E-01 

0.3045E-01 

0.2527E-01 

0.2646E-01 

0.2253E-01 

0.2330E-01 

0.1 81 2E— 01 

0.1838E-01 

0. 1667E-01 

0.1702E-01 

0.1173E-01 


0.9806E-02 


0.8838E-03 



0. 1000E+01 

0.5492E+00 

0.5492E+00 

0.4022E+00 

0.3559E+00 

0.3517E+00 

0.2999E+00 

0.2999E+00 

0.2627E+00 

0.2627E+00 

0.2501E+00 

0.2403E+00 

0.2394E+00 

0.2211E+00 

0.2209E+00 

0.2146E+00 

0.2146E+00 

0.2090E+00 

0.2089E+00 

0. 1943E+00 

0.1941E+00 

0. 1940E+00 

0. 1810E+00 

0. 1809E+00 

0. 1721E+00 


5.944E+04 


2.322E-02 


7.099E-02 


2.483E-02 


1 . 191 E - 02 
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TABLE 5.5 Summary of solution error for various lightning source geometries 
and locations. The effects of measurement error and constraints are tested. Note that 
solution charge errors are large, but that centroid errors are usually within the grid 
resolution ( = 2 km).* 
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Values of AX, AY, and Az represent spatial errors between the known centroid and the centroid of the distri- 
ition |f I. i.e., AX = ^retrieved * X true' ^ ~ ^ + ^ ^ s A * 2 + ^ 2 + A * 2 ; and AQ a (Z f - actual charge). 





























































TABLE 6.1 The total number of flashes identified during each storm period (N T ) , the 
number of flashes with Ae > 1 kV/m at 2 or more sites (N L ) , the number of 
acceptable solutions obtained with a Q— or P— model (N) , and the 
number of these solutions within the region of optimum accuracy (N Q ) . 
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